diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 6c4a5ce..78de918 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -12,7 +12,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        julia-version: ['1.6', '1', 'nightly']
+        julia-version: ['1', 'nightly']
         os: ['ubuntu-latest']
         include:
           - os: windows-latest
diff --git a/.github/workflows/docs-cleanup.yml b/.github/workflows/docs-cleanup.yml
index 0b00303..94b9f70 100644
--- a/.github/workflows/docs-cleanup.yml
+++ b/.github/workflows/docs-cleanup.yml
@@ -23,4 +23,4 @@ jobs:
               git push --force origin gh-pages-new:gh-pages
             fi
         env:
-            PRNUM: ${{ github.event.number }}
\ No newline at end of file
+            PRNUM: ${{ github.event.number }}
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index e423711..58db31b 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -9,12 +9,13 @@ on:
 jobs:
   build:
     runs-on: ubuntu-latest
+    timeout-minutes: 60
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@latest
         with:
-          version: '1.6'
-      - uses: actions/cache@v3
+          version: '1.10'
+      - uses: actions/cache@v4
         with:
           path: |
             ~/.julia/artifacts
@@ -22,7 +23,7 @@ jobs:
             ~/.julia/registries
           key: .julia-docs-${{ hashFiles('docs/Project.toml', 'docs/Manifest.toml') }}
       - name: Install dependencies
-        run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+        run: julia --project=docs -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()'
       - name: Build and deploy
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token
diff --git a/Project.toml b/Project.toml
index dd8eb59..997958d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,14 +8,16 @@ Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 ShaderAbstractions = "65257c39-d410-5151-9873-9b3e5be5013e"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 
 [compat]
-Ferrite = "0.3.13"
+Ferrite = "1"
 GeometryBasics = "0.4"
-Makie = "0.17,0.18,0.19"
-ShaderAbstractions = "0.3, 0.4"
+Makie = "0.21"
+ShaderAbstractions = "0.3,0.4"
 Tensors = "1"
 julia = "1"
 
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 2673e0a..7dc4fc1 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.9.3"
+julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "d59e06be2f3a2a869077caa7bc752d83cf9f5b47"
+project_hash = "3584c0ea6a470f2745c28be2dfe2f35e2c289a51"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
@@ -20,26 +20,32 @@ weakdeps = ["ChainRulesCore", "Test"]
     AbstractFFTsChainRulesCoreExt = "ChainRulesCore"
     AbstractFFTsTestExt = "Test"
 
-[[deps.AbstractLattices]]
-git-tree-sha1 = "f35684b7349da49fcc8a9e520e30e45dbb077166"
-uuid = "398f06c4-4d28-53ec-89ca-5b2656b7603d"
-version = "0.2.1"
-
 [[deps.AbstractTrees]]
-git-tree-sha1 = "faa260e4cb5aba097a73fab382dd4b5819d8ec8c"
+git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177"
 uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
-version = "0.4.4"
+version = "0.4.5"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "02f731463748db57cc2ebfbd9fbc9ce8280d3433"
+git-tree-sha1 = "50c3c56a52972d78e8be9fd135bfb91c9574c140"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.7.1"
+version = "4.1.1"
 weakdeps = ["StaticArrays"]
 
     [deps.Adapt.extensions]
     AdaptStaticArraysExt = "StaticArrays"
 
+[[deps.AdaptivePredicates]]
+git-tree-sha1 = "7e651ea8d262d2d74ce75fdf47c4d63c07dba7a6"
+uuid = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"
+version = "1.2.0"
+
+[[deps.AliasTables]]
+deps = ["PtrArrays", "Random"]
+git-tree-sha1 = "9876e1e164b144ca45e9e3198d0b689cadfed9ff"
+uuid = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"
+version = "1.1.3"
+
 [[deps.Animations]]
 deps = ["Colors"]
 git-tree-sha1 = "e81c509d2c8e49592413bfb0bb3b08150056c79d"
@@ -50,33 +56,11 @@ version = "0.4.1"
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 version = "1.1.1"
 
-[[deps.ArrayInterface]]
-deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "16267cf279190ca7c1b30d020758ced95db89cd0"
-uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.5.1"
-
-    [deps.ArrayInterface.extensions]
-    ArrayInterfaceBandedMatricesExt = "BandedMatrices"
-    ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
-    ArrayInterfaceCUDAExt = "CUDA"
-    ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
-    ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
-    ArrayInterfaceTrackerExt = "Tracker"
-
-    [deps.ArrayInterface.weakdeps]
-    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
-    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
-    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-    GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
-    StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
-
 [[deps.ArrayLayouts]]
 deps = ["FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "af43df5704827c8618afd36eb56fcab20d3041ee"
+git-tree-sha1 = "492681bc44fac86804706ddb37da10880a2bd528"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
-version = "1.4.3"
+version = "1.10.4"
 weakdeps = ["SparseArrays"]
 
     [deps.ArrayLayouts.extensions]
@@ -86,16 +70,16 @@ weakdeps = ["SparseArrays"]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
 [[deps.Automa]]
-deps = ["TranscodingStreams"]
-git-tree-sha1 = "ef9997b3d5547c48b41c7bd8899e812a917b409d"
+deps = ["PrecompileTools", "SIMD", "TranscodingStreams"]
+git-tree-sha1 = "a8f503e8e1a5f583fbef15a8440c8c7e32185df2"
 uuid = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b"
-version = "0.8.4"
+version = "1.1.0"
 
 [[deps.AxisAlgorithms]]
 deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
-git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7"
+git-tree-sha1 = "01b8ccb13d68535d73d2b0c23e39bd23155fb712"
 uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
-version = "1.0.1"
+version = "1.1.0"
 
 [[deps.AxisArrays]]
 deps = ["Dates", "IntervalSets", "IterTools", "RangeArrays"]
@@ -107,36 +91,42 @@ version = "0.4.7"
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
 [[deps.BitFlags]]
-git-tree-sha1 = "43b1a4a8f797c1cddadf60499a8a077d4af2cd2d"
+git-tree-sha1 = "0691e34b3bb8be9307330f88d1a3c3f25466c24d"
 uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35"
-version = "0.1.7"
+version = "0.1.9"
 
 [[deps.BlockArrays]]
 deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "54cd829dd26330c42e1cf9df68470dd4df602c61"
+git-tree-sha1 = "d434647f798823bcae510aee0bc0401927f64391"
 uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-version = "0.16.38"
+version = "1.1.1"
+
+    [deps.BlockArrays.extensions]
+    BlockArraysBandedMatricesExt = "BandedMatrices"
+
+    [deps.BlockArrays.weakdeps]
+    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+
+[[deps.Bonito]]
+deps = ["Base64", "CodecZlib", "Colors", "Dates", "Deno_jll", "HTTP", "Hyperscript", "LinearAlgebra", "Markdown", "MsgPack", "Observables", "RelocatableFolders", "SHA", "Sockets", "Tables", "ThreadPools", "URIs", "UUIDs", "WidgetsBase"]
+git-tree-sha1 = "d7635780a8cfe0cb43c075276fd358c5b166695e"
+uuid = "824d6782-a2ef-11e9-3a09-e5662e0c26f8"
+version = "3.2.4"
 
 [[deps.Bzip2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
+git-tree-sha1 = "8873e196c2eb87962a2048b3b8e08946535864a1"
 uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
-version = "1.0.8+0"
+version = "1.0.8+2"
 
 [[deps.CEnum]]
-git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90"
+git-tree-sha1 = "389ad5c84de1ae7cf0e28e381131c98ea87d54fc"
 uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
-version = "0.4.2"
+version = "0.5.0"
 
 [[deps.CRC32c]]
 uuid = "8bf52ea8-c179-5cab-976a-9e18b702a9bc"
 
-[[deps.CRlibm]]
-deps = ["CRlibm_jll"]
-git-tree-sha1 = "32abd86e3c2025db5172aa182b982debed519834"
-uuid = "96374032-68de-5a5b-8d9e-752f78720389"
-version = "1.0.1"
-
 [[deps.CRlibm_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "e329286945d0cfc04456972ea732551869af1cfc"
@@ -144,22 +134,16 @@ uuid = "4e9b3aee-d8a1-5a3d-ad8b-7d824db253f0"
 version = "1.0.1+0"
 
 [[deps.Cairo_jll]]
-deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
+deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
+git-tree-sha1 = "009060c9a6168704143100f36ab08f06c2af4642"
 uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
-version = "1.16.1+1"
-
-[[deps.Calculus]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
-uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
-version = "0.5.1"
+version = "1.18.2+1"
 
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra"]
-git-tree-sha1 = "e0af648f0692ec1691b5d094b8724ba1346281cf"
+git-tree-sha1 = "3e4b134270b372f2ed4d4d0e936aabaefc1802bc"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.18.0"
+version = "1.25.0"
 weakdeps = ["SparseArrays"]
 
     [deps.ChainRulesCore.extensions]
@@ -167,9 +151,9 @@ weakdeps = ["SparseArrays"]
 
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
-git-tree-sha1 = "cd67fc487743b2f0fd4380d4cbd3a24660d0eec8"
+git-tree-sha1 = "bce6804e5e6044c6daab27bb533d1295e4a2e759"
 uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
-version = "0.7.3"
+version = "0.7.6"
 
 [[deps.ColorBrewer]]
 deps = ["Colors", "JSON", "Test"]
@@ -179,44 +163,43 @@ version = "0.4.0"
 
 [[deps.ColorSchemes]]
 deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "PrecompileTools", "Random"]
-git-tree-sha1 = "67c1f244b991cad9b0aa4b7540fb758c2488b129"
+git-tree-sha1 = "13951eb68769ad1cd460cdb2e64e5e95f1bf123d"
 uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
-version = "3.24.0"
+version = "3.27.0"
 
 [[deps.ColorTypes]]
 deps = ["FixedPointNumbers", "Random"]
-git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4"
+git-tree-sha1 = "b10d0b65641d57b8b4d5e234446582de5047050d"
 uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
-version = "0.11.4"
+version = "0.11.5"
 
 [[deps.ColorVectorSpace]]
-deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
-git-tree-sha1 = "600cc5508d66b78aae350f7accdb58763ac18589"
+deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "Requires", "Statistics", "TensorCore"]
+git-tree-sha1 = "a1f44953f2382ebb937d60dafbe2deea4bd23249"
 uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
-version = "0.9.10"
+version = "0.10.0"
+weakdeps = ["SpecialFunctions"]
+
+    [deps.ColorVectorSpace.extensions]
+    SpecialFunctionsExt = "SpecialFunctions"
 
 [[deps.Colors]]
 deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
-git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a"
+git-tree-sha1 = "362a287c3aa50601b0bc359053d5c2468f0e7ce0"
 uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
-version = "0.12.10"
-
-[[deps.Combinatorics]]
-git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
-uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
-version = "1.0.2"
+version = "0.12.11"
 
 [[deps.CommonSubexpressions]]
-deps = ["MacroTools", "Test"]
-git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
+deps = ["MacroTools"]
+git-tree-sha1 = "cda2cfaebb4be89c9084adaca7dd7333369715c5"
 uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
-version = "0.3.0"
+version = "0.3.1"
 
 [[deps.Compat]]
-deps = ["UUIDs"]
-git-tree-sha1 = "8a62af3e248a8c4bad6b32cbbe663ae02275e32c"
+deps = ["TOML", "UUIDs"]
+git-tree-sha1 = "8ae8d32e09f0dcf42a36b90d4e17f5dd2e4c4215"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.10.0"
+version = "4.16.0"
 weakdeps = ["Dates", "LinearAlgebra"]
 
     [deps.Compat.extensions]
@@ -225,40 +208,40 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+0"
+version = "1.1.1+0"
 
 [[deps.ConcurrentUtilities]]
 deps = ["Serialization", "Sockets"]
-git-tree-sha1 = "8cfa272e8bdedfa88b6aefbbca7c19f1befac519"
+git-tree-sha1 = "ea32b83ca4fefa1768dc84e504cc0a94fb1ab8d1"
 uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb"
-version = "2.3.0"
+version = "2.4.2"
 
 [[deps.ConstructionBase]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "c53fc348ca4d40d7b371e71fd52251839080cbc9"
+git-tree-sha1 = "76219f1ed5771adbb096743bff43fb5fdd4c1157"
 uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
-version = "1.5.4"
-weakdeps = ["IntervalSets", "StaticArrays"]
+version = "1.5.8"
+weakdeps = ["IntervalSets", "LinearAlgebra", "StaticArrays"]
 
     [deps.ConstructionBase.extensions]
     ConstructionBaseIntervalSetsExt = "IntervalSets"
+    ConstructionBaseLinearAlgebraExt = "LinearAlgebra"
     ConstructionBaseStaticArraysExt = "StaticArrays"
 
 [[deps.Contour]]
-git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781"
+git-tree-sha1 = "439e35b0b36e2e5881738abc8857bd92ad6ff9a8"
 uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
-version = "0.6.2"
+version = "0.6.3"
 
 [[deps.DataAPI]]
-git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c"
+git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe"
 uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
-version = "1.15.0"
+version = "1.16.0"
 
 [[deps.DataStructures]]
 deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
-git-tree-sha1 = "3dbd312d370723b6bb43ba9d02fc36abade4518d"
+git-tree-sha1 = "1d0a14036acb104d9e89698bd408f63ab58cdc82"
 uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.18.15"
+version = "0.18.20"
 
 [[deps.DataValueInterfaces]]
 git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
@@ -270,10 +253,10 @@ deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
 [[deps.DelaunayTriangulation]]
-deps = ["DataStructures", "EnumX", "ExactPredicates", "Random", "SimpleGraphs"]
-git-tree-sha1 = "7cb0d72a53c1d93665eeadfa9d51af9df60bf6b2"
+deps = ["AdaptivePredicates", "EnumX", "ExactPredicates", "PrecompileTools", "Random"]
+git-tree-sha1 = "89df54fbe66e5872d91d8c2cd3a375f660c3fd64"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
-version = "0.8.10"
+version = "1.6.1"
 
 [[deps.Deno_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -295,9 +278,9 @@ version = "1.15.1"
 
 [[deps.Distances]]
 deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
-git-tree-sha1 = "5225c965635d8c21168e32a12954675e7bea1151"
+git-tree-sha1 = "c7e3a542b999843086e2f29dac96a618c105be1d"
 uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
-version = "0.10.10"
+version = "0.10.12"
 weakdeps = ["ChainRulesCore", "SparseArrays"]
 
     [deps.Distances.extensions]
@@ -309,10 +292,10 @@ deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
 [[deps.Distributions]]
-deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"]
-git-tree-sha1 = "a6c00f894f24460379cb7136633cef54ac9f6f4a"
+deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"]
+git-tree-sha1 = "d7477ecdafb813ddee2ae727afa94e9dcb5f3fb0"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.25.103"
+version = "0.25.112"
 
     [deps.Distributions.extensions]
     DistributionsChainRulesCoreExt = "ChainRulesCore"
@@ -331,22 +314,16 @@ uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 version = "0.9.3"
 
 [[deps.Documenter]]
-deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "Dates", "DocStringExtensions", "Downloads", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "Test", "Unicode"]
-git-tree-sha1 = "662fb21ae7fad33e044c2b59ece832fdce32c171"
+deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"]
+git-tree-sha1 = "5a1ee886566f2fa9318df1273d8b778b9d42712d"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "1.1.2"
+version = "1.7.0"
 
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
 uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 version = "1.6.0"
 
-[[deps.DualNumbers]]
-deps = ["Calculus", "NaNMath", "SpecialFunctions"]
-git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566"
-uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
-version = "0.6.8"
-
 [[deps.EarCut_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "e3290f2d49e661fbd94046d7e3726ffcb2d41053"
@@ -358,63 +335,52 @@ git-tree-sha1 = "bdb1942cd4c45e3c678fd11569d5cccd80976237"
 uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.4"
 
-[[deps.ErrorfreeArithmetic]]
-git-tree-sha1 = "d6863c556f1142a061532e79f611aa46be201686"
-uuid = "90fa49ef-747e-5e6f-a989-263ba693cf1a"
-version = "0.5.2"
-
 [[deps.ExactPredicates]]
-deps = ["IntervalArithmetic", "Random", "StaticArraysCore"]
-git-tree-sha1 = "499b1ca78f6180c8f8bdf1cabde2d39120229e5c"
+deps = ["IntervalArithmetic", "Random", "StaticArrays"]
+git-tree-sha1 = "b3f2ff58735b5f024c392fde763f29b057e4b025"
 uuid = "429591f6-91af-11e9-00e2-59fbe8cec110"
-version = "2.2.6"
+version = "2.2.8"
 
 [[deps.ExceptionUnwrapping]]
 deps = ["Test"]
-git-tree-sha1 = "e90caa41f5a86296e014e148ee061bd6c3edec96"
+git-tree-sha1 = "dcb08a0d93ec0b1cdc4af184b26b591e9695423a"
 uuid = "460bff9d-24e4-43bc-9d9f-a8973cb893f4"
-version = "0.1.9"
+version = "0.1.10"
 
 [[deps.Expat_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "4558ab818dcceaab612d1bb8c19cee87eda2b83c"
+git-tree-sha1 = "1c6317308b9dc757616f0b5cb379db10494443a7"
 uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
-version = "2.5.0+0"
+version = "2.6.2+0"
 
 [[deps.Extents]]
-git-tree-sha1 = "2140cd04483da90b2da7f99b2add0750504fc39c"
+git-tree-sha1 = "81023caa0021a41712685887db1fc03db26f41f5"
 uuid = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
-version = "0.1.2"
+version = "0.1.4"
 
 [[deps.FFMPEG_jll]]
 deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
-git-tree-sha1 = "466d45dc38e15794ec7d5d63ec03d776a9aff36e"
+git-tree-sha1 = "8cc47f299902e13f90405ddb5bf87e5d474c0d38"
 uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
-version = "4.4.4+1"
+version = "6.1.2+0"
 
 [[deps.FFTW]]
 deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
-git-tree-sha1 = "b4fbdd20c889804969571cc589900803edda16b7"
+git-tree-sha1 = "4820348781ae578893311153d69049a93d05f39d"
 uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-version = "1.7.1"
+version = "1.8.0"
 
 [[deps.FFTW_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea"
+git-tree-sha1 = "4d81ed14783ec49ce9f2e168208a12ce1815aa25"
 uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
-version = "3.3.10+0"
-
-[[deps.FastRounding]]
-deps = ["ErrorfreeArithmetic", "LinearAlgebra"]
-git-tree-sha1 = "6344aa18f654196be82e62816935225b3b9abe44"
-uuid = "fa42c844-2597-5d31-933b-ebd51ab2693f"
-version = "0.3.1"
+version = "3.3.10+1"
 
 [[deps.Ferrite]]
-deps = ["EnumX", "LinearAlgebra", "NearestNeighbors", "Preferences", "Reexport", "SparseArrays", "Tensors", "WriteVTK"]
-git-tree-sha1 = "a31b9d4dd58e00686e5da175dab11667af3108b6"
+deps = ["EnumX", "ForwardDiff", "LinearAlgebra", "NearestNeighbors", "OrderedCollections", "Preferences", "Reexport", "SparseArrays", "StaticArrays", "Tensors", "WriteVTK"]
+git-tree-sha1 = "1ad6e7eab1803998ad16360b1f86c7723d361cb0"
 uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
-version = "0.3.14"
+version = "1.0.0"
 
     [deps.Ferrite.extensions]
     FerriteBlockArrays = "BlockArrays"
@@ -425,70 +391,71 @@ version = "0.3.14"
     Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
 
 [[deps.FerriteViz]]
-deps = ["Ferrite", "GeometryBasics", "LinearAlgebra", "Makie", "ShaderAbstractions", "Tensors"]
+deps = ["Ferrite", "GeometryBasics", "LinearAlgebra", "Makie", "Reexport", "ShaderAbstractions", "StaticArrays", "Tensors"]
 path = ".."
 uuid = "59d0093e-b1f1-4fb7-ac85-ab57e45f39d9"
-version = "0.2.1"
+version = "0.2.2"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
-git-tree-sha1 = "299dc33549f68299137e51e6d49a13b5b1da9673"
+git-tree-sha1 = "62ca0547a14c57e98154423419d8a342dca75ca9"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.16.1"
+version = "1.16.4"
+
+[[deps.FilePaths]]
+deps = ["FilePathsBase", "MacroTools", "Reexport", "Requires"]
+git-tree-sha1 = "919d9412dbf53a2e6fe74af62a73ceed0bce0629"
+uuid = "8fc22ac5-c921-52a6-82fd-178b2807b824"
+version = "0.8.3"
+
+[[deps.FilePathsBase]]
+deps = ["Compat", "Dates"]
+git-tree-sha1 = "7878ff7172a8e6beedd1dea14bd27c3c6340d361"
+uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
+version = "0.9.22"
+weakdeps = ["Mmap", "Test"]
+
+    [deps.FilePathsBase.extensions]
+    FilePathsBaseMmapExt = "Mmap"
+    FilePathsBaseTestExt = "Test"
 
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
 [[deps.FillArrays]]
-deps = ["LinearAlgebra", "Random"]
-git-tree-sha1 = "35f0c0f345bff2c6d636f95fdb136323b5a796ef"
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "6a70198746448456524cb442b8af316927ff3e1a"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.7.0"
-weakdeps = ["SparseArrays", "Statistics"]
+version = "1.13.0"
+weakdeps = ["PDMats", "SparseArrays", "Statistics"]
 
     [deps.FillArrays.extensions]
+    FillArraysPDMatsExt = "PDMats"
     FillArraysSparseArraysExt = "SparseArrays"
     FillArraysStatisticsExt = "Statistics"
 
-[[deps.FiniteDiff]]
-deps = ["ArrayInterface", "LinearAlgebra", "Requires", "Setfield", "SparseArrays"]
-git-tree-sha1 = "c6e4a1fbe73b31a3dea94b1da449503b8830c306"
-uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.21.1"
-
-    [deps.FiniteDiff.extensions]
-    FiniteDiffBandedMatricesExt = "BandedMatrices"
-    FiniteDiffBlockBandedMatricesExt = "BlockBandedMatrices"
-    FiniteDiffStaticArraysExt = "StaticArrays"
-
-    [deps.FiniteDiff.weakdeps]
-    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
-    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
-    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-
 [[deps.FixedPointNumbers]]
 deps = ["Statistics"]
-git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
+git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172"
 uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
-version = "0.8.4"
+version = "0.8.5"
 
 [[deps.Fontconfig_jll]]
-deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
+deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"]
+git-tree-sha1 = "db16beca600632c95fc8aca29890d83788dd8b23"
 uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
-version = "2.13.93+0"
+version = "2.13.96+0"
 
-[[deps.Formatting]]
-deps = ["Printf"]
-git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
-uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
-version = "0.4.2"
+[[deps.Format]]
+git-tree-sha1 = "9c68794ef81b08086aeb32eeaf33531668d5f5fc"
+uuid = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
+version = "1.3.7"
 
 [[deps.ForwardDiff]]
 deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"]
-git-tree-sha1 = "cf0fe81336da9fb90944683b8c41984b08793dad"
+git-tree-sha1 = "a9ce73d3c827adab2d70bf168aaece8cce196898"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.36"
+version = "0.10.37"
 weakdeps = ["StaticArrays"]
 
     [deps.ForwardDiff.extensions]
@@ -496,49 +463,44 @@ weakdeps = ["StaticArrays"]
 
 [[deps.FreeType]]
 deps = ["CEnum", "FreeType2_jll"]
-git-tree-sha1 = "50351f83f95282cf903e968d7c6e8d44a5f83d0b"
+git-tree-sha1 = "907369da0f8e80728ab49c1c7e09327bf0d6d999"
 uuid = "b38be410-82b0-50bf-ab77-7b57e271db43"
-version = "4.1.0"
+version = "4.1.1"
 
 [[deps.FreeType2_jll]]
 deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "d8db6a5a2fe1381c1ea4ef2cab7c69c2de7f9ea0"
+git-tree-sha1 = "5c1d8ae0efc6c2e7b1fc502cbe25def8f661b7bc"
 uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
-version = "2.13.1+0"
+version = "2.13.2+0"
 
 [[deps.FreeTypeAbstraction]]
 deps = ["ColorVectorSpace", "Colors", "FreeType", "GeometryBasics"]
-git-tree-sha1 = "38a92e40157100e796690421e34a11c107205c86"
+git-tree-sha1 = "84dfe824bd6fdf2a5d73bb187ff31b5549b2a79c"
 uuid = "663a7486-cb36-511b-a19d-713bb74d65c9"
-version = "0.10.0"
+version = "0.10.4"
 
 [[deps.FriBidi_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "1ed150b39aebcc805c26b93a8d0122c940f64ce2"
 uuid = "559328eb-81f9-559d-9380-de523a88c83c"
-version = "1.0.10+0"
+version = "1.0.14+0"
 
-[[deps.Future]]
-deps = ["Random"]
-uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
-
-[[deps.GPUArraysCore]]
-deps = ["Adapt"]
-git-tree-sha1 = "2d6ca471a6c7b536127afccfa7564b5b39227fe0"
-uuid = "46192b85-c4d5-4398-a991-12ede77f4527"
-version = "0.1.5"
+[[deps.GeoFormatTypes]]
+git-tree-sha1 = "59107c179a586f0fe667024c5eb7033e81333271"
+uuid = "68eda718-8dee-11e9-39e7-89f7f65f511f"
+version = "0.4.2"
 
 [[deps.GeoInterface]]
-deps = ["Extents"]
-git-tree-sha1 = "d53480c0793b13341c40199190f92c611aa2e93c"
+deps = ["Extents", "GeoFormatTypes"]
+git-tree-sha1 = "2f6fce56cdb8373637a6614e14a5768a88450de2"
 uuid = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
-version = "1.3.2"
+version = "1.3.7"
 
 [[deps.GeometryBasics]]
 deps = ["EarCut_jll", "Extents", "GeoInterface", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
-git-tree-sha1 = "424a5a6ce7c5d97cca7bcc4eac551b97294c54af"
+git-tree-sha1 = "b62f2b2d76cee0d61a2ef2b3118cd2a3215d3134"
 uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
-version = "0.4.9"
+version = "0.4.11"
 
 [[deps.Gettext_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
@@ -546,17 +508,29 @@ git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
 uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
 version = "0.21.0+0"
 
+[[deps.Giflib_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "0224cce99284d997f6880a42ef715a37c99338d1"
+uuid = "59f7168a-df46-5410-90c8-f2779963d0ec"
+version = "5.2.2+0"
+
+[[deps.Git]]
+deps = ["Git_jll"]
+git-tree-sha1 = "04eff47b1354d702c3a85e8ab23d539bb7d5957e"
+uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2"
+version = "1.3.1"
+
+[[deps.Git_jll]]
+deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"]
+git-tree-sha1 = "ea372033d09e4552a04fd38361cd019f9003f4f4"
+uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
+version = "2.46.2+0"
+
 [[deps.Glib_jll]]
 deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "e94c92c7bf4819685eb80186d51c43e71d4afa17"
+git-tree-sha1 = "674ff0db93fffcd11a3573986e550d66cd4fd71f"
 uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
-version = "2.76.5+0"
-
-[[deps.Graphics]]
-deps = ["Colors", "LinearAlgebra", "NaNMath"]
-git-tree-sha1 = "d61890399bc535850c4bf08e4e0d3a7ad0f21cbd"
-uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
-version = "1.1.2"
+version = "2.80.5+0"
 
 [[deps.Graphite2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -566,9 +540,9 @@ version = "1.3.14+0"
 
 [[deps.GridLayoutBase]]
 deps = ["GeometryBasics", "InteractiveUtils", "Observables"]
-git-tree-sha1 = "f57a64794b336d4990d90f80b147474b869b1bc4"
+git-tree-sha1 = "fc713f007cff99ff9e50accba6373624ddd33588"
 uuid = "3955a311-db13-416c-9275-1d80ed98e5e9"
-version = "0.9.2"
+version = "0.11.0"
 
 [[deps.Grisu]]
 git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
@@ -577,33 +551,33 @@ version = "1.0.2"
 
 [[deps.HTTP]]
 deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"]
-git-tree-sha1 = "5eab648309e2e060198b45820af1a37182de3cce"
+git-tree-sha1 = "bc3f416a965ae61968c20d0ad867556367f2817d"
 uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
-version = "1.10.0"
+version = "1.10.9"
 
 [[deps.HarfBuzz_jll]]
-deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
-git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
+deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll"]
+git-tree-sha1 = "401e4f3f30f43af2c8478fc008da50096ea5240f"
 uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
-version = "2.8.1+1"
+version = "8.3.1+0"
 
 [[deps.HypergeometricFunctions]]
-deps = ["DualNumbers", "LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"]
-git-tree-sha1 = "f218fe3736ddf977e0e772bc9a586b2383da2685"
+deps = ["LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"]
+git-tree-sha1 = "7c4195be1649ae622304031ed46a2f4df989f1eb"
 uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
-version = "0.3.23"
+version = "0.3.24"
 
 [[deps.Hyperscript]]
 deps = ["Test"]
-git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
+git-tree-sha1 = "179267cfa5e712760cd43dcae385d7ea90cc25a4"
 uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
-version = "0.0.4"
+version = "0.0.5"
 
 [[deps.IOCapture]]
 deps = ["Logging", "Random"]
-git-tree-sha1 = "d75853a0bdbfb1ac815478bacd89cd27b550ace6"
+git-tree-sha1 = "b6d6bfdd7ce25b0f9b2f6b3dd56b2673a66c8770"
 uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
-version = "0.2.3"
+version = "0.2.5"
 
 [[deps.ImageAxes]]
 deps = ["AxisArrays", "ImageBase", "ImageCore", "Reexport", "SimpleTraits"]
@@ -613,21 +587,21 @@ version = "0.6.11"
 
 [[deps.ImageBase]]
 deps = ["ImageCore", "Reexport"]
-git-tree-sha1 = "b51bb8cae22c66d0f6357e3bcb6363145ef20835"
+git-tree-sha1 = "eb49b82c172811fd2c86759fa0553a2221feb909"
 uuid = "c817782e-172a-44cc-b673-b171935fbb9e"
-version = "0.1.5"
+version = "0.1.7"
 
 [[deps.ImageCore]]
-deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"]
-git-tree-sha1 = "acf614720ef026d38400b3817614c45882d75500"
+deps = ["ColorVectorSpace", "Colors", "FixedPointNumbers", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "PrecompileTools", "Reexport"]
+git-tree-sha1 = "b2a7eaa169c13f5bcae8131a83bc30eff8f71be0"
 uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
-version = "0.9.4"
+version = "0.10.2"
 
 [[deps.ImageIO]]
-deps = ["FileIO", "IndirectArrays", "JpegTurbo", "LazyModules", "Netpbm", "OpenEXR", "PNGFiles", "QOI", "Sixel", "TiffImages", "UUIDs"]
-git-tree-sha1 = "bca20b2f5d00c4fbc192c3212da8fa79f4688009"
+deps = ["FileIO", "IndirectArrays", "JpegTurbo", "LazyModules", "Netpbm", "OpenEXR", "PNGFiles", "QOI", "Sixel", "TiffImages", "UUIDs", "WebP"]
+git-tree-sha1 = "696144904b76e1ca433b886b4e7edd067d76cbf7"
 uuid = "82e4d734-157c-48bb-816b-45c225c6df19"
-version = "0.6.7"
+version = "0.6.9"
 
 [[deps.ImageMetadata]]
 deps = ["AxisArrays", "ImageAxes", "ImageBase", "ImageCore"]
@@ -637,9 +611,9 @@ version = "0.9.9"
 
 [[deps.Imath_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "3d09a9f60edf77f8a4d99f9e015e8fbf9989605d"
+git-tree-sha1 = "0936ba688c6d201805a83da835b55c61a180db52"
 uuid = "905a6f67-0a94-5f89-b386-d35d92009cd1"
-version = "3.1.7+0"
+version = "3.1.11+0"
 
 [[deps.IndirectArrays]]
 git-tree-sha1 = "012e604e1c7458645cb8b436f8fba789a51b257f"
@@ -647,20 +621,15 @@ uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
 version = "1.0.0"
 
 [[deps.Inflate]]
-git-tree-sha1 = "ea8031dea4aff6bd41f1df8f2fdfb25b33626381"
+git-tree-sha1 = "d1b1b796e47d94588b3757fe84fbf65a5ec4a80d"
 uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
-version = "0.1.4"
-
-[[deps.IntegerMathUtils]]
-git-tree-sha1 = "b8ffb903da9f7b8cf695a8bead8e01814aa24b30"
-uuid = "18e54dd8-cb9d-406c-a71d-865a43cbb235"
-version = "0.1.2"
+version = "0.1.5"
 
 [[deps.IntelOpenMP_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "ad37c091f7d7daf900963171600d7c1c5c3ede32"
+deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"]
+git-tree-sha1 = "10bd689145d2c3b2a9844005d01087cc1194e79e"
 uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
-version = "2023.2.0+0"
+version = "2024.2.1+0"
 
 [[deps.InteractiveUtils]]
 deps = ["Markdown"]
@@ -668,26 +637,48 @@ uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [[deps.Interpolations]]
 deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
-git-tree-sha1 = "721ec2cf720536ad005cb38f50dbba7b02419a15"
+git-tree-sha1 = "88a101217d7cb38a7b481ccd50d21876e1d1b0e0"
 uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
-version = "0.14.7"
+version = "0.15.1"
+weakdeps = ["Unitful"]
+
+    [deps.Interpolations.extensions]
+    InterpolationsUnitfulExt = "Unitful"
 
 [[deps.IntervalArithmetic]]
-deps = ["CRlibm", "EnumX", "FastRounding", "LinearAlgebra", "Markdown", "Random", "RecipesBase", "RoundingEmulator", "SetRounding", "StaticArrays"]
-git-tree-sha1 = "f59e639916283c1d2e106d2b00910b50f4dab76c"
+deps = ["CRlibm_jll", "LinearAlgebra", "MacroTools", "RoundingEmulator"]
+git-tree-sha1 = "c59c57c36683aa17c563be6edaac888163f35285"
 uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
-version = "0.21.2"
+version = "0.22.18"
+weakdeps = ["DiffRules", "ForwardDiff", "IntervalSets", "RecipesBase"]
+
+    [deps.IntervalArithmetic.extensions]
+    IntervalArithmeticDiffRulesExt = "DiffRules"
+    IntervalArithmeticForwardDiffExt = "ForwardDiff"
+    IntervalArithmeticIntervalSetsExt = "IntervalSets"
+    IntervalArithmeticRecipesBaseExt = "RecipesBase"
 
 [[deps.IntervalSets]]
-deps = ["Dates", "Random"]
-git-tree-sha1 = "3d8866c029dd6b16e69e0d4a939c4dfcb98fac47"
+git-tree-sha1 = "dba9ddf07f77f60450fe5d2e2beb9854d9a49bd0"
 uuid = "8197267c-284f-5f27-9208-e0e47529a953"
-version = "0.7.8"
-weakdeps = ["Statistics"]
+version = "0.7.10"
+weakdeps = ["Random", "RecipesBase", "Statistics"]
 
     [deps.IntervalSets.extensions]
+    IntervalSetsRandomExt = "Random"
+    IntervalSetsRecipesBaseExt = "RecipesBase"
     IntervalSetsStatisticsExt = "Statistics"
 
+[[deps.InverseFunctions]]
+git-tree-sha1 = "a779299d77cd080bf77b97535acecd73e1c5e5cb"
+uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
+version = "0.1.17"
+weakdeps = ["Dates", "Test"]
+
+    [deps.InverseFunctions.extensions]
+    InverseFunctionsDatesExt = "Dates"
+    InverseFunctionsTestExt = "Test"
+
 [[deps.IrrationalConstants]]
 git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
 uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
@@ -700,9 +691,9 @@ uuid = "f1662d9f-8043-43de-a69a-05efc1cc6ff4"
 version = "0.1.1"
 
 [[deps.IterTools]]
-git-tree-sha1 = "4ced6667f9974fc5c5943fa5e2ef1ca43ea9e450"
+git-tree-sha1 = "42d5f897009e7ff2cf88db414a389e5ed1bdd023"
 uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
-version = "1.8.0"
+version = "1.10.0"
 
 [[deps.IteratorInterfaceExtensions]]
 git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
@@ -711,9 +702,9 @@ version = "1.0.0"
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca"
+git-tree-sha1 = "be3dc50a92e5a386872a493a10050136d4703f9b"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.5.0"
+version = "1.6.1"
 
 [[deps.JSON]]
 deps = ["Dates", "Mmap", "Parsers", "Unicode"]
@@ -721,57 +712,57 @@ git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 version = "0.21.4"
 
-[[deps.JSServe]]
-deps = ["Base64", "CodecZlib", "Colors", "Dates", "Deno_jll", "HTTP", "Hyperscript", "LinearAlgebra", "Markdown", "MsgPack", "Observables", "RelocatableFolders", "SHA", "Sockets", "Tables", "ThreadPools", "URIs", "UUIDs", "WidgetsBase"]
-git-tree-sha1 = "2f849fa6a48753c5c715e9d630405fd2d0855573"
-uuid = "824d6782-a2ef-11e9-3a09-e5662e0c26f9"
-version = "2.3.0"
-
 [[deps.JpegTurbo]]
 deps = ["CEnum", "FileIO", "ImageCore", "JpegTurbo_jll", "TOML"]
-git-tree-sha1 = "d65930fa2bc96b07d7691c652d701dcbe7d9cf0b"
+git-tree-sha1 = "fa6d0bcff8583bac20f1ffa708c3913ca605c611"
 uuid = "b835a17e-a41a-41e7-81f0-2f016b05efe0"
-version = "0.1.4"
+version = "0.1.5"
 
 [[deps.JpegTurbo_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "6f2675ef130a300a112286de91973805fcc5ffbc"
+git-tree-sha1 = "25ee0be4d43d0269027024d75a24c24d6c6e590c"
 uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
-version = "2.1.91+0"
+version = "3.0.4+0"
 
 [[deps.KernelDensity]]
 deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"]
-git-tree-sha1 = "90442c50e202a5cdf21a7899c66b240fdef14035"
+git-tree-sha1 = "7d703202e65efa1369de1279c162b915e245eed1"
 uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
-version = "0.6.7"
+version = "0.6.9"
 
 [[deps.LAME_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "170b660facf5df5de098d866564877e119141cbd"
 uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
-version = "3.100.1+0"
+version = "3.100.2+0"
+
+[[deps.LERC_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "36bdbc52f13a7d1dcb0f3cd694e01677a515655b"
+uuid = "88015f11-f218-50d7-93a8-a6af411a945d"
+version = "4.0.0+0"
 
 [[deps.LLVMOpenMP_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "f689897ccbe049adb19a065c495e75f372ecd42b"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "78211fb6cbc872f77cad3fc0b6cf647d923f4929"
 uuid = "1d63c593-3942-5779-bab2-d838dc0a180e"
-version = "15.0.4+0"
+version = "18.1.7+0"
 
 [[deps.LZO_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "854a9c268c43b77b0a27f22d7fab8d33cdb3a731"
 uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
-version = "2.10.1+0"
+version = "2.10.2+1"
 
 [[deps.LaTeXStrings]]
-git-tree-sha1 = "50901ebc375ed41dbf8058da26f9de442febbbec"
+git-tree-sha1 = "dda21b8cbd6a6c40d9d02a73230f9d70fed6918c"
 uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
-version = "1.3.1"
+version = "1.4.0"
 
 [[deps.LazilyInitializedFields]]
-git-tree-sha1 = "410fe4739a4b092f2ffe36fcb0dcc3ab12648ce1"
+git-tree-sha1 = "0f2da712350b020bc3957f269c9caad516383ee0"
 uuid = "0e77f7df-68c5-4e49-93ce-4cd80f5598bf"
-version = "1.2.1"
+version = "1.3.0"
 
 [[deps.LazyArtifacts]]
 deps = ["Artifacts", "Pkg"]
@@ -785,21 +776,26 @@ version = "0.3.1"
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
-version = "0.6.3"
+version = "0.6.4"
 
 [[deps.LibCURL_jll]]
 deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
 uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
-version = "7.84.0+0"
+version = "8.4.0+0"
 
 [[deps.LibGit2]]
-deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
+[[deps.LibGit2_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"]
+uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5"
+version = "1.6.4+0"
+
 [[deps.LibSSH2_jll]]
 deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
 uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
-version = "1.10.2+0"
+version = "1.11.0+1"
 
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -811,68 +807,68 @@ uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
 version = "3.2.2+1"
 
 [[deps.Libgcrypt_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
-git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll"]
+git-tree-sha1 = "8be878062e0ffa2c3f67bb58a595375eda5de80b"
 uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
-version = "1.8.7+0"
+version = "1.11.0+0"
+
+[[deps.Libglvnd_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
+git-tree-sha1 = "6f73d1dd803986947b2c750138528a999a6c7733"
+uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
+version = "1.6.0+0"
 
 [[deps.Libgpg_error_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "c6ce1e19f3aec9b59186bdf06cdf3c4fc5f5f3e6"
 uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
-version = "1.42.0+0"
+version = "1.50.0+0"
 
 [[deps.Libiconv_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175"
+git-tree-sha1 = "61dfdba58e585066d8bce214c5a51eaa0539f269"
 uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
-version = "1.17.0+0"
+version = "1.17.0+1"
 
 [[deps.Libmount_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "0c4f9c4f1a50d8f35048fa0532dabbadf702f81e"
 uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
-version = "2.35.0+0"
+version = "2.40.1+0"
+
+[[deps.Libtiff_jll]]
+deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
+git-tree-sha1 = "b404131d06f7886402758c9ce2214b636eb4d54a"
+uuid = "89763e89-9b03-5906-acba-b20f662cd828"
+version = "4.7.0+0"
 
 [[deps.Libuuid_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "5ee6203157c120d79034c748a2acba45b82b8807"
 uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
-version = "2.36.0+0"
+version = "2.40.1+0"
 
 [[deps.LightXML]]
 deps = ["Libdl", "XML2_jll"]
-git-tree-sha1 = "e129d9391168c677cd4800f5c0abb1ed8cb3794f"
+git-tree-sha1 = "3a994404d3f6709610701c7dabfc03fed87a81f8"
 uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
-version = "0.9.0"
-
-[[deps.LineSearches]]
-deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
-git-tree-sha1 = "7bbea35cec17305fc70a0e5b4641477dc0789d9d"
-uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
-version = "7.2.0"
+version = "0.9.1"
 
 [[deps.LinearAlgebra]]
 deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
-[[deps.LinearAlgebraX]]
-deps = ["LinearAlgebra", "Mods", "Permutations", "Primes", "SimplePolynomials"]
-git-tree-sha1 = "558a338f1eeabe933f9c2d4052aa7c2c707c3d52"
-uuid = "9b3f67b0-2d00-526e-9884-9e4938f8fb88"
-version = "0.1.12"
-
 [[deps.LiveServer]]
-deps = ["HTTP", "LoggingExtras", "MIMEs", "Pkg", "Sockets", "Test"]
-git-tree-sha1 = "24d05efe53436b22a42bf2ae459f47c48b0c2603"
+deps = ["HTTP", "LoggingExtras", "MIMEs", "Sockets", "Test"]
+git-tree-sha1 = "564a436267fb1fc768f815dad64c4386c46623f8"
 uuid = "16fef848-5104-11e9-1b77-fb7a48bbb589"
-version = "1.2.7"
+version = "1.4.0"
 
 [[deps.LogExpFunctions]]
 deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa"
+git-tree-sha1 = "a2d09619db4e765091ee5c6ffe8872849de0feea"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.26"
+version = "0.3.28"
 
     [deps.LogExpFunctions.extensions]
     LogExpFunctionsChainRulesCoreExt = "ChainRulesCore"
@@ -889,9 +885,9 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
 [[deps.LoggingExtras]]
 deps = ["Dates", "Logging"]
-git-tree-sha1 = "c1dd6d7978c12545b4179fb6153b9250c96b0075"
+git-tree-sha1 = "f02b56007b064fbfddb4c9cd60161b6dd0f40df3"
 uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36"
-version = "1.0.3"
+version = "1.1.0"
 
 [[deps.MIMEs]]
 git-tree-sha1 = "65f28ad4b594aebe22157d6fac869786a255b7eb"
@@ -899,28 +895,28 @@ uuid = "6c6e2e6c-3030-632d-7369-2d6c69616d65"
 version = "0.1.4"
 
 [[deps.MKL_jll]]
-deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
-git-tree-sha1 = "eb006abbd7041c28e0d16260e50a24f8f9104913"
+deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"]
+git-tree-sha1 = "f046ccd0c6db2832a9f639e2c669c6fe867e5f4f"
 uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
-version = "2023.2.0+0"
+version = "2024.2.0+0"
 
 [[deps.MacroTools]]
 deps = ["Markdown", "Random"]
-git-tree-sha1 = "9ee1618cbf5240e6d4e0371d6f24065083f60c48"
+git-tree-sha1 = "2fa9ee3e63fd3a4f7a9a4f4744a52f4856de82df"
 uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-version = "0.5.11"
+version = "0.5.13"
 
 [[deps.Makie]]
-deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FixedPointNumbers", "Formatting", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Setfield", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "StableHashTraits", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun"]
-git-tree-sha1 = "35fa3c150cd96fd77417a23965b7037b90d6ffc9"
+deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "Dates", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageBase", "ImageIO", "InteractiveUtils", "Interpolations", "IntervalSets", "InverseFunctions", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun", "Unitful"]
+git-tree-sha1 = "f7907907eb914138cc9e9ee66ab46f7a9efac8e8"
 uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
-version = "0.19.12"
+version = "0.21.15"
 
 [[deps.MakieCore]]
-deps = ["Observables", "REPL"]
-git-tree-sha1 = "9b11acd07f21c4d035bd4156e789532e8ee2cc70"
+deps = ["ColorTypes", "GeometryBasics", "IntervalSets", "Observables"]
+git-tree-sha1 = "4604f03e5b057e8e62a95a44929cafc9585b0fe9"
 uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
-version = "0.6.9"
+version = "0.8.9"
 
 [[deps.MappedArrays]]
 git-tree-sha1 = "2dab0221fe2b0f2cb6754eaa743cc266339f527e"
@@ -938,36 +934,31 @@ uuid = "d0879d2d-cac2-40c8-9cee-1863dc0c7391"
 version = "0.1.2"
 
 [[deps.MathTeXEngine]]
-deps = ["AbstractTrees", "Automa", "DataStructures", "FreeTypeAbstraction", "GeometryBasics", "LaTeXStrings", "REPL", "RelocatableFolders", "Test", "UnicodeFun"]
-git-tree-sha1 = "8f52dbaa1351ce4cb847d95568cb29e62a307d93"
+deps = ["AbstractTrees", "Automa", "DataStructures", "FreeTypeAbstraction", "GeometryBasics", "LaTeXStrings", "REPL", "RelocatableFolders", "UnicodeFun"]
+git-tree-sha1 = "f45c8916e8385976e1ccd055c9874560c257ab13"
 uuid = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
-version = "0.5.6"
+version = "0.6.2"
 
 [[deps.MbedTLS]]
-deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "Random", "Sockets"]
-git-tree-sha1 = "03a9b9718f5682ecb107ac9f7308991db4ce395b"
+deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "NetworkOptions", "Random", "Sockets"]
+git-tree-sha1 = "c067a280ddc25f196b5e7df3877c6b226d390aaf"
 uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
-version = "1.1.7"
+version = "1.1.9"
 
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
-version = "2.28.2+0"
+version = "2.28.2+1"
 
 [[deps.Missings]]
 deps = ["DataAPI"]
-git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272"
+git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d"
 uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
-version = "1.1.0"
+version = "1.2.0"
 
 [[deps.Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
-[[deps.Mods]]
-git-tree-sha1 = "61be59e4daffff43a8cec04b5e0dc773cbb5db3a"
-uuid = "7475f97c-0381-53b1-977b-4c60186c8d62"
-version = "1.3.3"
-
 [[deps.MosaicViews]]
 deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"]
 git-tree-sha1 = "7b86a5d4d70a9f5cdf2dacb3cbe6d251d1a61dbe"
@@ -976,24 +967,13 @@ version = "0.3.4"
 
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
-version = "2022.10.11"
+version = "2023.1.10"
 
 [[deps.MsgPack]]
 deps = ["Serialization"]
-git-tree-sha1 = "fc8c15ca848b902015bd4a745d350f02cf791c2a"
+git-tree-sha1 = "f5db02ae992c260e4826fe78c942954b48e1d9c2"
 uuid = "99f44e22-a591-53d1-9472-aa23ef4bd671"
-version = "1.2.0"
-
-[[deps.Multisets]]
-git-tree-sha1 = "8d852646862c96e226367ad10c8af56099b4047e"
-uuid = "3b2b4ff1-bcff-5658-a3ee-dbcf1ce5ac09"
-version = "0.4.4"
-
-[[deps.NLSolversBase]]
-deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
-git-tree-sha1 = "a0b464d183da839699f4c79e7606d9d186ec172c"
-uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
-version = "7.8.3"
+version = "1.2.1"
 
 [[deps.NaNMath]]
 deps = ["OpenLibm_jll"]
@@ -1003,9 +983,9 @@ version = "1.0.2"
 
 [[deps.NearestNeighbors]]
 deps = ["Distances", "StaticArrays"]
-git-tree-sha1 = "2c3726ceb3388917602169bed973dbc97f1b51a8"
+git-tree-sha1 = "3cebfc94a0754cc329ebc3bab1e6c89621e791ad"
 uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
-version = "0.4.13"
+version = "0.4.20"
 
 [[deps.Netpbm]]
 deps = ["FileIO", "ImageCore", "ImageMetadata"]
@@ -1023,10 +1003,13 @@ uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
 version = "0.5.5"
 
 [[deps.OffsetArrays]]
-deps = ["Adapt"]
-git-tree-sha1 = "2ac17d29c523ce1cd38e27785a7d23024853a4bb"
+git-tree-sha1 = "1a27764e945a152f7ca7efa04de513d473e9542e"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.12.10"
+version = "1.14.1"
+weakdeps = ["Adapt"]
+
+    [deps.OffsetArrays.extensions]
+    OffsetArraysAdaptExt = "Adapt"
 
 [[deps.Ogg_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -1037,7 +1020,7 @@ version = "1.3.5+1"
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.21+4"
+version = "0.3.23+4"
 
 [[deps.OpenEXR]]
 deps = ["Colors", "FileIO", "OpenEXR_jll"]
@@ -1047,26 +1030,26 @@ version = "0.3.2"
 
 [[deps.OpenEXR_jll]]
 deps = ["Artifacts", "Imath_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "a4ca623df1ae99d09bc9868b008262d0c0ac1e4f"
+git-tree-sha1 = "8292dd5c8a38257111ada2174000a33745b06d4e"
 uuid = "18a262bb-aa17-5467-a713-aee519bc75cb"
-version = "3.1.4+0"
+version = "3.2.4+0"
 
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.1+0"
+version = "0.8.1+2"
 
 [[deps.OpenSSL]]
 deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"]
-git-tree-sha1 = "51901a49222b09e3743c65b8847687ae5fc78eb2"
+git-tree-sha1 = "38cb508d080d21dc1128f7fb04f20387ed4c0af4"
 uuid = "4d8831e6-92b7-49fb-bdf8-b643e874388c"
-version = "1.4.1"
+version = "1.4.3"
 
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "cc6e1927ac521b659af340e0ca45828a3ffc748f"
+git-tree-sha1 = "7493f61f55a6cce7325f197443aa80d32554ba10"
 uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
-version = "3.0.12+0"
+version = "3.0.15+1"
 
 [[deps.OpenSpecFun_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
@@ -1074,39 +1057,33 @@ git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
 uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 version = "0.5.5+0"
 
-[[deps.Optim]]
-deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"]
-git-tree-sha1 = "01f85d9269b13fedc61e63cc72ee2213565f7a72"
-uuid = "429524aa-4258-5aef-a3af-852621145aeb"
-version = "1.7.8"
-
 [[deps.Opus_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "6703a85cb3781bd5909d48730a67205f3f31a575"
 uuid = "91d4177d-7536-5919-b921-800302f37372"
-version = "1.3.2+0"
+version = "1.3.3+0"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "2e73fe17cac3c62ad1aebe70d44c963c3cfdc3e3"
+git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.6.2"
+version = "1.6.3"
 
 [[deps.PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15"
-version = "10.42.0+0"
+version = "10.42.0+1"
 
 [[deps.PDMats]]
 deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "66b2fcd977db5329aa35cac121e5b94dd6472198"
+git-tree-sha1 = "949347156c25054de2db3b166c52ac4728cbad65"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.11.28"
+version = "0.11.31"
 
 [[deps.PNGFiles]]
 deps = ["Base64", "CEnum", "ImageCore", "IndirectArrays", "OffsetArrays", "libpng_jll"]
-git-tree-sha1 = "5ded86ccaf0647349231ed6c0822c10886d4a1ee"
+git-tree-sha1 = "67186a2bc9a90f9f85ff3cc8277868961fb57cbd"
 uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
-version = "0.4.1"
+version = "0.4.3"
 
 [[deps.Packing]]
 deps = ["GeometryBasics"]
@@ -1120,34 +1097,22 @@ git-tree-sha1 = "0fac6313486baae819364c52b4f483450a9d793f"
 uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
 version = "0.5.12"
 
-[[deps.Parameters]]
-deps = ["OrderedCollections", "UnPack"]
-git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
-uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
-version = "0.12.3"
-
 [[deps.Parsers]]
 deps = ["Dates", "PrecompileTools", "UUIDs"]
-git-tree-sha1 = "716e24b21538abc91f6205fd1d8363f39b442851"
+git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.7.2"
-
-[[deps.Permutations]]
-deps = ["Combinatorics", "LinearAlgebra", "Random"]
-git-tree-sha1 = "4f69b02cf40a0f494d0438ab29de32e14ef96e7b"
-uuid = "2ae35dd2-176d-5d53-8349-f30d82d94d4f"
-version = "0.4.18"
+version = "2.8.1"
 
 [[deps.Pixman_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"]
-git-tree-sha1 = "64779bc4c9784fee475689a1752ef4d5747c5e87"
+git-tree-sha1 = "35621f10a7531bc8fa58f74610b1bfb70a3cfc6b"
 uuid = "30392449-352a-5448-841d-b1acce4e97dc"
-version = "0.42.2+0"
+version = "0.43.4+0"
 
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-version = "1.9.2"
+version = "1.10.0"
 
 [[deps.PkgVersion]]
 deps = ["Pkg"]
@@ -1156,57 +1121,27 @@ uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
 version = "0.3.3"
 
 [[deps.PlotUtils]]
-deps = ["ColorSchemes", "Colors", "Dates", "PrecompileTools", "Printf", "Random", "Reexport", "Statistics"]
-git-tree-sha1 = "f92e1315dadf8c46561fb9396e525f7200cdc227"
+deps = ["ColorSchemes", "Colors", "Dates", "PrecompileTools", "Printf", "Random", "Reexport", "StableRNGs", "Statistics"]
+git-tree-sha1 = "3ca9a356cd2e113c420f2c13bea19f8d3fb1cb18"
 uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
-version = "1.3.5"
+version = "1.4.3"
 
 [[deps.PolygonOps]]
 git-tree-sha1 = "77b3d3605fc1cd0b42d95eba87dfcd2bf67d5ff6"
 uuid = "647866c9-e3ac-4575-94e7-e3d426903924"
 version = "0.1.2"
 
-[[deps.Polynomials]]
-deps = ["LinearAlgebra", "RecipesBase", "Setfield", "SparseArrays"]
-git-tree-sha1 = "ea78a2764f31715093de7ab495e12c0187f231d1"
-uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
-version = "4.0.4"
-
-    [deps.Polynomials.extensions]
-    PolynomialsChainRulesCoreExt = "ChainRulesCore"
-    PolynomialsFFTWExt = "FFTW"
-    PolynomialsMakieCoreExt = "MakieCore"
-    PolynomialsMutableArithmeticsExt = "MutableArithmetics"
-
-    [deps.Polynomials.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-    FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-    MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
-    MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
-
-[[deps.PositiveFactorizations]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20"
-uuid = "85a6dd25-e78a-55b7-8502-1745935b8125"
-version = "0.2.4"
-
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
+git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.2.0"
+version = "1.2.1"
 
 [[deps.Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e"
+git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.4.1"
-
-[[deps.Primes]]
-deps = ["IntegerMathUtils"]
-git-tree-sha1 = "4c9f306e5d6603ae203c2000dd460d81a5251489"
-uuid = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
-version = "0.5.4"
+version = "1.4.3"
 
 [[deps.Printf]]
 deps = ["Unicode"]
@@ -1214,9 +1149,14 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [[deps.ProgressMeter]]
 deps = ["Distributed", "Printf"]
-git-tree-sha1 = "00099623ffee15972c16111bcf84c58a0051257c"
+git-tree-sha1 = "8f6bc219586aef8baf0ff9a5fe16ee9c70cb65e4"
 uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
-version = "1.9.0"
+version = "1.10.2"
+
+[[deps.PtrArrays]]
+git-tree-sha1 = "77a42d78b6a92df47ab37e177b2deac405e1c88f"
+uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d"
+version = "1.2.1"
 
 [[deps.QOI]]
 deps = ["ColorTypes", "FileIO", "FixedPointNumbers"]
@@ -1226,16 +1166,22 @@ version = "1.0.0"
 
 [[deps.QuadGK]]
 deps = ["DataStructures", "LinearAlgebra"]
-git-tree-sha1 = "9ebcd48c498668c7fa0e97a9cae873fbee7bfee1"
+git-tree-sha1 = "cda3b045cf9ef07a08ad46731f5a3165e56cf3da"
 uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
-version = "2.9.1"
+version = "2.11.1"
+
+    [deps.QuadGK.extensions]
+    QuadGKEnzymeExt = "Enzyme"
+
+    [deps.QuadGK.weakdeps]
+    Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 
 [[deps.REPL]]
 deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 
 [[deps.Random]]
-deps = ["SHA", "Serialization"]
+deps = ["SHA"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[deps.RangeArrays]]
@@ -1282,23 +1228,17 @@ git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
 version = "1.3.0"
 
-[[deps.RingLists]]
-deps = ["Random"]
-git-tree-sha1 = "f39da63aa6d2d88e0c1bd20ed6a3ff9ea7171ada"
-uuid = "286e9d63-9694-5540-9e3c-4e6708fa07b2"
-version = "0.2.8"
-
 [[deps.Rmath]]
 deps = ["Random", "Rmath_jll"]
-git-tree-sha1 = "f65dcb5fa46aee0cf9ed6274ccbd597adc49aa7b"
+git-tree-sha1 = "852bd0f55565a9e973fcfee83a84413270224dc4"
 uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
-version = "0.7.1"
+version = "0.8.0"
 
 [[deps.Rmath_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "6ed52fdd3382cf21947b15e8870ac0ddbff736da"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "58cdd8fb2201a6267e1db87ff148dd6c1dbd8ad8"
 uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
-version = "0.4.0+0"
+version = "0.5.1+0"
 
 [[deps.RoundingEmulator]]
 git-tree-sha1 = "40b9edad2e5287e05bd413a38f61a8ff55b9557b"
@@ -1311,9 +1251,9 @@ version = "0.7.0"
 
 [[deps.SIMD]]
 deps = ["PrecompileTools"]
-git-tree-sha1 = "d8911cc125da009051fb35322415641d02d9e37f"
+git-tree-sha1 = "98ca7c29edd6fc79cd74c61accb7010a4e7aee33"
 uuid = "fdea26ae-647d-5447-a871-4b548cad5224"
-version = "3.4.6"
+version = "3.6.0"
 
 [[deps.Scratch]]
 deps = ["Dates"]
@@ -1324,22 +1264,11 @@ version = "1.2.1"
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
-[[deps.SetRounding]]
-git-tree-sha1 = "d7a25e439d07a17b7cdf97eecee504c50fedf5f6"
-uuid = "3cc68bcd-71a2-5612-b932-767ffbe40ab0"
-version = "0.2.1"
-
-[[deps.Setfield]]
-deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"]
-git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac"
-uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
-version = "1.1.1"
-
 [[deps.ShaderAbstractions]]
 deps = ["ColorTypes", "FixedPointNumbers", "GeometryBasics", "LinearAlgebra", "Observables", "StaticArrays", "StructArrays", "Tables"]
-git-tree-sha1 = "db0219befe4507878b1a90e07820fed3e62c289d"
+git-tree-sha1 = "79123bc60c5507f035e6d1d9e563bb2971954ec8"
 uuid = "65257c39-d410-5151-9873-9b3e5be5013e"
-version = "0.4.0"
+version = "0.4.1"
 
 [[deps.SharedArrays]]
 deps = ["Distributed", "Mmap", "Random", "Serialization"]
@@ -1358,33 +1287,9 @@ uuid = "73760f76-fbc4-59ce-8f25-708e95d2df96"
 version = "0.4.0"
 
 [[deps.SimpleBufferStream]]
-git-tree-sha1 = "874e8867b33a00e784c8a7e4b60afe9e037b74e1"
+git-tree-sha1 = "f305871d2f381d21527c770d4788c06c097c9bc1"
 uuid = "777ac1f9-54b0-4bf8-805c-2214025038e7"
-version = "1.1.0"
-
-[[deps.SimpleGraphs]]
-deps = ["AbstractLattices", "Combinatorics", "DataStructures", "IterTools", "LightXML", "LinearAlgebra", "LinearAlgebraX", "Optim", "Primes", "Random", "RingLists", "SimplePartitions", "SimplePolynomials", "SimpleRandom", "SparseArrays", "Statistics"]
-git-tree-sha1 = "b608903049d11cc557c45e03b3a53e9260579c19"
-uuid = "55797a34-41de-5266-9ec1-32ac4eb504d3"
-version = "0.8.4"
-
-[[deps.SimplePartitions]]
-deps = ["AbstractLattices", "DataStructures", "Permutations"]
-git-tree-sha1 = "dcc02923a53f316ab97da8ef3136e80b4543dbf1"
-uuid = "ec83eff0-a5b5-5643-ae32-5cbf6eedec9d"
-version = "0.3.0"
-
-[[deps.SimplePolynomials]]
-deps = ["Mods", "Multisets", "Polynomials", "Primes"]
-git-tree-sha1 = "d537c31cf9995236166e3e9afc424a5a1c59ff9d"
-uuid = "cc47b68c-3164-5771-a705-2bc0097375a0"
-version = "0.2.14"
-
-[[deps.SimpleRandom]]
-deps = ["Distributions", "LinearAlgebra", "Random"]
-git-tree-sha1 = "3a6fb395e37afab81aeea85bae48a4db5cd7244a"
-uuid = "a6525b86-64cd-54fa-8f65-62fc48bdc0e8"
-version = "0.3.1"
+version = "1.2.0"
 
 [[deps.SimpleTraits]]
 deps = ["InteractiveUtils", "MacroTools"]
@@ -1403,29 +1308,30 @@ uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
 [[deps.SortingAlgorithms]]
 deps = ["DataStructures"]
-git-tree-sha1 = "5165dfb9fd131cf0c6957a3a7605dede376e7b63"
+git-tree-sha1 = "66e0a8e672a0bdfca2c3f5937efb8538b9ddc085"
 uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
-version = "1.2.0"
+version = "1.2.1"
 
 [[deps.SparseArrays]]
 deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+version = "1.10.0"
 
 [[deps.SpecialFunctions]]
 deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d"
+git-tree-sha1 = "2f5d4697f21388cbe1ff299430dd169ef97d7e14"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.3.1"
+version = "2.4.0"
 weakdeps = ["ChainRulesCore"]
 
     [deps.SpecialFunctions.extensions]
     SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
 
-[[deps.StableHashTraits]]
-deps = ["Compat", "SHA", "Tables", "TupleTools"]
-git-tree-sha1 = "d29023a76780bb8a3f2273b29153fd00828cb73f"
-uuid = "c5dd0088-6c3f-4803-b00e-f31a60c170fa"
-version = "1.1.1"
+[[deps.StableRNGs]]
+deps = ["Random"]
+git-tree-sha1 = "83e6cce8324d49dfaf9ef059227f91ed4441a8e5"
+uuid = "860ef19b-820b-49d6-a774-d7a799459cd3"
+version = "1.0.2"
 
 [[deps.StackViews]]
 deps = ["OffsetArrays"]
@@ -1434,24 +1340,25 @@ uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
 version = "0.1.1"
 
 [[deps.StaticArrays]]
-deps = ["LinearAlgebra", "Random", "StaticArraysCore"]
-git-tree-sha1 = "0adf069a2a490c47273727e029371b31d44b72b2"
+deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
+git-tree-sha1 = "777657803913ffc7e8cc20f0fd04b634f871af8f"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.6.5"
-weakdeps = ["Statistics"]
+version = "1.9.8"
+weakdeps = ["ChainRulesCore", "Statistics"]
 
     [deps.StaticArrays.extensions]
+    StaticArraysChainRulesCoreExt = "ChainRulesCore"
     StaticArraysStatisticsExt = "Statistics"
 
 [[deps.StaticArraysCore]]
-git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
+git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682"
 uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-version = "1.4.2"
+version = "1.4.3"
 
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-version = "1.9.0"
+version = "1.10.0"
 
 [[deps.StatsAPI]]
 deps = ["LinearAlgebra"]
@@ -1461,38 +1368,47 @@ version = "1.7.0"
 
 [[deps.StatsBase]]
 deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
-git-tree-sha1 = "1d77abd07f617c4868c33d4f5b9e1dbb2643c9cf"
+git-tree-sha1 = "5cf7606d6cef84b543b483848d4ae08ad9832b21"
 uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-version = "0.34.2"
+version = "0.34.3"
 
 [[deps.StatsFuns]]
 deps = ["HypergeometricFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
-git-tree-sha1 = "f625d686d5a88bcd2b15cd81f18f98186fdc0c9a"
+git-tree-sha1 = "b423576adc27097764a90e163157bcfc9acf0f46"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "1.3.0"
+version = "1.3.2"
+weakdeps = ["ChainRulesCore", "InverseFunctions"]
 
     [deps.StatsFuns.extensions]
     StatsFunsChainRulesCoreExt = "ChainRulesCore"
     StatsFunsInverseFunctionsExt = "InverseFunctions"
 
-    [deps.StatsFuns.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-    InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
-
 [[deps.StructArrays]]
-deps = ["Adapt", "ConstructionBase", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"]
-git-tree-sha1 = "0a3db38e4cce3c54fe7a71f831cd7b6194a54213"
+deps = ["ConstructionBase", "DataAPI", "Tables"]
+git-tree-sha1 = "f4dc295e983502292c4c3f951dbb4e985e35b3be"
 uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
-version = "0.6.16"
+version = "0.6.18"
+
+    [deps.StructArrays.extensions]
+    StructArraysAdaptExt = "Adapt"
+    StructArraysGPUArraysCoreExt = "GPUArraysCore"
+    StructArraysSparseArraysExt = "SparseArrays"
+    StructArraysStaticArraysExt = "StaticArrays"
+
+    [deps.StructArrays.weakdeps]
+    Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+    GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [[deps.SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
 uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [[deps.SuiteSparse_jll]]
-deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"]
+deps = ["Artifacts", "Libdl", "libblastrampoline_jll"]
 uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
-version = "5.10.1+6"
+version = "7.2.1+1"
 
 [[deps.TOML]]
 deps = ["Dates"]
@@ -1506,10 +1422,10 @@ uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
 version = "1.0.1"
 
 [[deps.Tables]]
-deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits"]
-git-tree-sha1 = "cb76cf677714c095e535e3501ac7954732aeea2d"
+deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "OrderedCollections", "TableTraits"]
+git-tree-sha1 = "598cd7c1f68d1e205689b1c2fe65a9f85846f297"
 uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "1.11.1"
+version = "1.12.0"
 
 [[deps.Tar]]
 deps = ["ArgTools", "SHA"]
@@ -1539,27 +1455,21 @@ uuid = "b189fb0b-2eb5-4ed4-bc0c-d34c51242431"
 version = "2.1.1"
 
 [[deps.TiffImages]]
-deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "Mmap", "OffsetArrays", "PkgVersion", "ProgressMeter", "UUIDs"]
-git-tree-sha1 = "34cc045dd0aaa59b8bbe86c644679bc57f1d5bd0"
+deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "Mmap", "OffsetArrays", "PkgVersion", "ProgressMeter", "SIMD", "UUIDs"]
+git-tree-sha1 = "6ee0c220d0aecad18792c277ae358129cc50a475"
 uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69"
-version = "0.6.8"
+version = "0.11.0"
 
 [[deps.TranscodingStreams]]
-deps = ["Random", "Test"]
-git-tree-sha1 = "9a6ae7ed916312b41236fcef7e0af564ef934769"
+git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
-version = "0.9.13"
+version = "0.11.3"
 
 [[deps.TriplotBase]]
 git-tree-sha1 = "4d4ed7f294cda19382ff7de4c137d24d16adc89b"
 uuid = "981d1d27-644d-49a2-9326-4793e63143c3"
 version = "0.1.0"
 
-[[deps.TupleTools]]
-git-tree-sha1 = "155515ed4c4236db30049ac1495e2969cc06be9d"
-uuid = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
-version = "1.4.3"
-
 [[deps.URIs]]
 git-tree-sha1 = "67db6cc7b3821e19ebe75791a9dd19c9b1188f2b"
 uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
@@ -1569,11 +1479,6 @@ version = "1.5.1"
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
-[[deps.UnPack]]
-git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
-uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
-version = "1.0.2"
-
 [[deps.Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
@@ -1583,16 +1488,33 @@ git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
 uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
 version = "0.4.1"
 
+[[deps.Unitful]]
+deps = ["Dates", "LinearAlgebra", "Random"]
+git-tree-sha1 = "d95fe458f26209c66a187b1114df96fd70839efd"
+uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
+version = "1.21.0"
+weakdeps = ["ConstructionBase", "InverseFunctions"]
+
+    [deps.Unitful.extensions]
+    ConstructionBaseUnitfulExt = "ConstructionBase"
+    InverseFunctionsUnitfulExt = "InverseFunctions"
+
 [[deps.VTKBase]]
 git-tree-sha1 = "c2d0db3ef09f1942d08ea455a9e252594be5f3b6"
 uuid = "4004b06d-e244-455f-a6ce-a5f9919cc534"
 version = "1.0.1"
 
 [[deps.WGLMakie]]
-deps = ["Colors", "FileIO", "FreeTypeAbstraction", "GeometryBasics", "Hyperscript", "JSServe", "LinearAlgebra", "Makie", "Observables", "PNGFiles", "PrecompileTools", "RelocatableFolders", "ShaderAbstractions", "StaticArrays"]
-git-tree-sha1 = "2328be9950072c7250d340798b25222779ce55bb"
+deps = ["Bonito", "Colors", "FileIO", "FreeTypeAbstraction", "GeometryBasics", "Hyperscript", "LinearAlgebra", "Makie", "Observables", "PNGFiles", "PrecompileTools", "RelocatableFolders", "ShaderAbstractions", "StaticArrays"]
+git-tree-sha1 = "13cab94d885d7760d487d7e2f30c2d6df4643880"
 uuid = "276b4fcb-3e11-5398-bf8b-a0c2d153d008"
-version = "0.8.16"
+version = "0.10.15"
+
+[[deps.WebP]]
+deps = ["CEnum", "ColorTypes", "FileIO", "FixedPointNumbers", "ImageCore", "libwebp_jll"]
+git-tree-sha1 = "f1f6d497ff84039deeb37f264396dac0c2250497"
+uuid = "e3aaa7dc-3e4b-44e0-be63-ffb868ccd7c1"
+version = "0.1.2"
 
 [[deps.WidgetsBase]]
 deps = ["Observables"]
@@ -1602,27 +1524,33 @@ version = "0.1.4"
 
 [[deps.WoodburyMatrices]]
 deps = ["LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "de67fa59e33ad156a590055375a30b23c40299d3"
+git-tree-sha1 = "c1a7aa6219628fcd757dede0ca95e245c5cd9511"
 uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
-version = "0.5.5"
+version = "1.0.0"
 
 [[deps.WriteVTK]]
 deps = ["Base64", "CodecZlib", "FillArrays", "LightXML", "TranscodingStreams", "VTKBase"]
-git-tree-sha1 = "41f0dc2a8f6fd860c266b91fd5cdf4fead65ae69"
+git-tree-sha1 = "1d8042d58334ab7947ce505709df7009da6f3375"
 uuid = "64499a7a-5c06-52f2-abe2-ccb03c286192"
-version = "1.18.1"
+version = "1.21.1"
 
 [[deps.XML2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
-git-tree-sha1 = "24b81b59bd35b3c42ab84fa589086e19be919916"
+git-tree-sha1 = "6a451c6f33a176150f315726eba8b92fbfdb9ae7"
 uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
-version = "2.11.5+0"
+version = "2.13.4+0"
 
 [[deps.XSLT_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
-git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "XML2_jll", "Zlib_jll"]
+git-tree-sha1 = "a54ee957f4c86b526460a720dbc882fa5edcbefc"
 uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
-version = "1.1.34+0"
+version = "1.1.41+0"
+
+[[deps.XZ_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "15e637a697345f6743674f1322beefbc5dcd5cfc"
+uuid = "ffd25f8a-64ca-5728-b0f7-c24cf3aae800"
+version = "5.6.3+0"
 
 [[deps.Xorg_libX11_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
@@ -1643,16 +1571,16 @@ uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
 version = "1.1.4+0"
 
 [[deps.Xorg_libXext_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
-git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
+git-tree-sha1 = "d2d1a5c49fae4ba39983f63de6afcbea47194e85"
 uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
-version = "1.3.4+4"
+version = "1.3.6+0"
 
 [[deps.Xorg_libXrender_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
-git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
+git-tree-sha1 = "47e45cd78224c53109495b3e324df0c37bb61fbe"
 uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
-version = "0.9.10+4"
+version = "0.9.11+0"
 
 [[deps.Xorg_libpthread_stubs_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1662,9 +1590,9 @@ version = "0.1.1+0"
 
 [[deps.Xorg_libxcb_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
-git-tree-sha1 = "b4bfde5d5b652e22b9c790ad00af08b6d042b97d"
+git-tree-sha1 = "bcd466676fef0878338c61e655629fa7bbc69d8e"
 uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
-version = "1.15.0+0"
+version = "1.17.0+0"
 
 [[deps.Xorg_xtrans_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1675,7 +1603,13 @@ version = "1.5.0+0"
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
-version = "1.2.13+0"
+version = "1.2.13+1"
+
+[[deps.Zstd_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "555d1076590a6cc2fdee2ef1469451f872d8b41b"
+uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
+version = "1.5.6+1"
 
 [[deps.isoband_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -1684,64 +1618,76 @@ uuid = "9a68df92-36a6-505f-a73e-abb412b6bfb4"
 version = "0.2.3+0"
 
 [[deps.libaom_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "3a2ea60308f0996d26f1e5354e10c24e9ef905d4"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "1827acba325fdcdf1d2647fc8d5301dd9ba43a9d"
 uuid = "a4ae2306-e953-59d6-aa16-d00cac43593b"
-version = "3.4.0+0"
+version = "3.9.0+0"
 
 [[deps.libass_jll]]
-deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
+deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
+git-tree-sha1 = "e17c115d55c5fbb7e52ebedb427a0dca79d4484e"
 uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
-version = "0.15.1+0"
+version = "0.15.2+0"
 
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.8.0+0"
+version = "5.8.0+1"
 
 [[deps.libfdk_aac_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "8a22cf860a7d27e4f3498a0fe0811a7957badb38"
 uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
-version = "2.0.2+0"
+version = "2.0.3+0"
 
 [[deps.libpng_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"]
+git-tree-sha1 = "b70c870239dc3d7bc094eb2d6be9b73d27bef280"
 uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
-version = "1.6.38+0"
+version = "1.6.44+0"
 
 [[deps.libsixel_jll]]
 deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "libpng_jll"]
-git-tree-sha1 = "d4f63314c8aa1e48cd22aa0c17ed76cd1ae48c3c"
+git-tree-sha1 = "7dfa0fd9c783d3d0cc43ea1af53d69ba45c447df"
 uuid = "075b6546-f08a-558a-be8f-8157d0f608a5"
-version = "1.10.3+0"
+version = "1.10.3+1"
 
 [[deps.libvorbis_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
-git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c"
+git-tree-sha1 = "490376214c4721cdaca654041f635213c6165cb3"
 uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
-version = "1.3.7+1"
+version = "1.3.7+2"
+
+[[deps.libwebp_jll]]
+deps = ["Artifacts", "Giflib_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libglvnd_jll", "Libtiff_jll", "libpng_jll"]
+git-tree-sha1 = "ccbb625a89ec6195856a50aa2b668a5c08712c94"
+uuid = "c5f90fcd-3b7e-5836-afba-fc50a0988cb2"
+version = "1.4.0+0"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
-version = "1.48.0+0"
+version = "1.52.0+1"
+
+[[deps.oneTBB_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "7d0ea0f4895ef2f5cb83645fa689e52cb55cf493"
+uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e"
+version = "2021.12.0+0"
 
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
-version = "17.4.0+0"
+version = "17.4.0+2"
 
 [[deps.x264_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "35976a1216d6c066ea32cba2150c4fa682b276fc"
 uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
-version = "2021.5.5+0"
+version = "10164.0.0+0"
 
 [[deps.x265_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "dcc541bb19ed5b0ede95581fb2e41ecf179527d2"
 uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
-version = "3.5.0+0"
+version = "3.6.0+0"
diff --git a/docs/Project.toml b/docs/Project.toml
index 926847d..93e7ce0 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,9 +1,9 @@
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+Bonito = "824d6782-a2ef-11e9-3a09-e5662e0c26f8"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
 FerriteViz = "59d0093e-b1f1-4fb7-ac85-ab57e45f39d9"
-JSServe = "824d6782-a2ef-11e9-3a09-e5662e0c26f9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
diff --git a/docs/make.jl b/docs/make.jl
index 1d3b445..47d13e7 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -2,6 +2,8 @@ push!(LOAD_PATH,"../src/")
 using Documenter, FerriteViz, WGLMakie
 
 makedocs(sitename="FerriteViz",
+         doctest = true,
+         warnonly = true,
          modules=[FerriteViz],
          authors="Maximilian Köhler",
          format=Documenter.HTML(
@@ -11,8 +13,8 @@ makedocs(sitename="FerriteViz",
          pages=["Home"=> "index.md",
                 "Tutorial" => "tutorial.md",
                 "Advanced Topics" => "atopics.md",
-                "api.md",
-                "devdocs.md",
+                "API Reference" => "api.md",
+                "Devdocs" => "devdocs.md",
                 ],
 )
 
diff --git a/docs/src/api.md b/docs/src/api.md
index 4bf0527..f452780 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -5,17 +5,11 @@ On this page the docs of the provided functions are listed
 ```@docs
 FerriteViz.MakiePlotter
 FerriteViz.solutionplot
-FerriteViz.solutionplot!
 FerriteViz.cellplot
-FerriteViz.cellplot!
 FerriteViz.wireframe
-FerriteViz.wireframe!
 FerriteViz.arrows
-FerriteViz.arrows!
 FerriteViz.surface
-FerriteViz.surface!
 FerriteViz.elementinfo
-FerriteViz.elementinfo!
 FerriteViz.ferriteviewer
 FerriteViz.update!
 FerriteViz.for_discretization
diff --git a/docs/src/atopics.md b/docs/src/atopics.md
index 99157fe..43c01ea 100644
--- a/docs/src/atopics.md
+++ b/docs/src/atopics.md
@@ -1,8 +1,10 @@
 # Advanced Topics
 
 ```@example 1
-import JSServe # hide
-JSServe.Page() # hide
+import WGLMakie, Bonito # hide
+Bonito.Page() # hide
+WGLMakie.activate!() # hide
+WGLMakie.Makie.inline!(true) # hide
 ```
 
 ## Gradient field visualization
@@ -17,7 +19,7 @@ In this quick example we show how to visualize strains and stresses side-by-side
 ```@example 1
 using Ferrite
 import FerriteViz
-using FerriteViz: ε
+ε(∇u) = (∇u+transpose(∇u))/2
 import WGLMakie #activating the backend, switch to GLMakie or CairoMakie (for 2D) locally
 
 include("ferrite-examples/incompressible-elasticity.jl") #defines dh_linear, dh_quadratic, u_linear, u_quadratic and mp
@@ -56,19 +58,19 @@ An alternative to this approach is to compute gradient quantities at samples poi
 The investigation of high-order fields is currently only supported via a first-order refinment of the problem.
 Here, the high-order approximation is replaced by a first order approximation of the field, which is
 spanned by the nodes of the high-order approximation. For example, the first order refinement of a
-heat problem on a square domain for Lagrange polynomials of order 5 looks like this:
+heat problem on a square domain for Lagrange polynomials of order 4 looks like this:
 ```@example 1
 include("ferrite-examples/heat-equation.jl"); #defines manufactured_heat_problem
 
 f = WGLMakie.Figure()
 axs = [WGLMakie.Axis3(f[1, 1], title="Coarse"), WGLMakie.Axis3(f[1, 2], title="Fine")]
 
-dh,u = manufactured_heat_problem(Triangle, Lagrange{2,RefTetrahedron,5}(), 1)
+dh,u = manufactured_heat_problem(Triangle, Lagrange{RefTriangle,4}(), 1)
 dh_for,u_for = FerriteViz.for_discretization(dh, u)
 plotter_for = FerriteViz.MakiePlotter(dh_for, u_for)
 FerriteViz.surface!(axs[1], plotter_for)
 
-dh,u = manufactured_heat_problem(Triangle, Lagrange{2,RefTetrahedron,5}(), 3)
+dh,u = manufactured_heat_problem(Triangle, Lagrange{RefTriangle,4}(), 3)
 dh_for,u_for = FerriteViz.for_discretization(dh, u)
 plotter_for = FerriteViz.MakiePlotter(dh_for, u_for)
 FerriteViz.surface!(axs[2], plotter_for)
@@ -84,7 +86,7 @@ include("ferrite-examples/heat-equation.jl"); #defines manufactured_heat_problem
 f = WGLMakie.Figure()
 axs = [WGLMakie.Axis3(f[1, 1], title="Coarse"), WGLMakie.Axis3(f[1, 2], title="Fine")]
 
-dh, u = manufactured_heat_problem(Hexahedron, Lagrange{3,RefCube,2}(), 2);
+dh, u = manufactured_heat_problem(Hexahedron, Lagrange{RefHexahedron,2}(), 2);
 plotter = FerriteViz.MakiePlotter(dh,u);
 clip_plane = FerriteViz.ClipPlane(Ferrite.Vec((0.0,0.5,0.5)), 0.1);
 clipped_plotter = FerriteViz.crinkle_clip(plotter, clip_plane);
diff --git a/docs/src/devdocs.md b/docs/src/devdocs.md
index e1c7e00..05be1b3 100644
--- a/docs/src/devdocs.md
+++ b/docs/src/devdocs.md
@@ -4,11 +4,9 @@ Note that these functions could be removed or change in behavior between minor v
 
 ```@docs
 FerriteViz.num_vertices
+FerriteViz.ntriangles
 FerriteViz.vertices
-FerriteViz.transfer_quadrature_face_to_cell
-FerriteViz.decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Union{FerriteViz.Ferrite.AbstractCell{2,N,3}, FerriteViz.Ferrite.AbstractCell{3,3,1}}) where {N}
-FerriteViz.decompose!(coord_offset, coord_matrix::Vector{Point{space_dim,T}}, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Union{FerriteViz.Ferrite.AbstractCell{2,N,4}, FerriteViz.Ferrite.AbstractCell{3,4,1}}) where {N,space_dim,T}
-FerriteViz.decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::FerriteViz.Ferrite.AbstractCell{3,N,M}) where {N,M}
+FerriteViz.decompose!
 FerriteViz.transfer_solution
 FerriteViz.postprocess
 FerriteViz._tensorsjl_gradient_accessor
diff --git a/docs/src/ferrite-examples/heat-equation.jl b/docs/src/ferrite-examples/heat-equation.jl
index 6acd069..1c5d036 100644
--- a/docs/src/ferrite-examples/heat-equation.jl
+++ b/docs/src/ferrite-examples/heat-equation.jl
@@ -1,6 +1,6 @@
 using Ferrite, SparseArrays
 
-function assemble_heat_element!(Ke::Matrix, fe::Vector, cellvalues::CellScalarValues, coords::Vector, rhs::Function)
+function assemble_heat_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues, coords::Vector, rhs::Function)
     n_basefuncs = getnbasefunctions(cellvalues)
 
     fill!(Ke, 0)
@@ -25,7 +25,7 @@ function assemble_heat_element!(Ke::Matrix, fe::Vector, cellvalues::CellScalarVa
     return Ke, fe
 end
 
-function assemble_steady_heat_global(cellvalues::CellScalarValues, K::SparseMatrixCSC, dh::DofHandler, rhs::Function)
+function assemble_steady_heat_global(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler, rhs::Function)
     n_basefuncs = getnbasefunctions(cellvalues)
     Ke = zeros(n_basefuncs, n_basefuncs)
     fe = zeros(n_basefuncs)
@@ -43,32 +43,32 @@ function assemble_steady_heat_global(cellvalues::CellScalarValues, K::SparseMatr
 end
 
 function manufactured_heat_problem(element_type, ip, num_elements_per_dim)
-    dim = Ferrite.getdim(ip)
+    dim = Ferrite.getrefdim(ip)
     grid = generate_grid(element_type, ntuple(x->num_elements_per_dim, dim));
-    ip_geo = Ferrite.default_interpolation(typeof(grid.cells[1]))
-    qr = QuadratureRule{dim, Ferrite.getrefshape(ip)}(2*Ferrite.getorder(ip))
-    cellvalues = CellScalarValues(qr, ip, ip_geo);
+    ip_geo = Ferrite.geometric_interpolation(typeof(grid.cells[1]))
+    qr = QuadratureRule{Ferrite.getrefshape(ip)}(2*Ferrite.getorder(ip))
+    cellvalues = CellValues(qr, ip, ip_geo);
 
     ∂Ω = union(
-        getfaceset(grid, "left"),
-        getfaceset(grid, "right"),
-        getfaceset(grid, "top"),
-        getfaceset(grid, "bottom"),
+        getfacetset(grid, "left"),
+        getfacetset(grid, "right"),
+        getfacetset(grid, "top"),
+        getfacetset(grid, "bottom"),
     );
 
     if dim == 3
         ∂Ω = union(
             ∂Ω,
-            getfaceset(grid, "front"),
-            getfaceset(grid, "back")
+            getfacetset(grid, "front"),
+            getfacetset(grid, "back")
         )
     end
 
     dh = DofHandler(grid)
-    push!(dh, :u, 1, ip)
+    add!(dh, :u, ip)
     close!(dh);
 
-    K = create_sparsity_pattern(dh)
+    K = allocate_matrix(dh)
 
     ch = ConstraintHandler(dh);
     dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0)
diff --git a/docs/src/ferrite-examples/incompressible-elasticity.jl b/docs/src/ferrite-examples/incompressible-elasticity.jl
index 6633176..5f1db49 100644
--- a/docs/src/ferrite-examples/incompressible-elasticity.jl
+++ b/docs/src/ferrite-examples/incompressible-elasticity.jl
@@ -8,40 +8,40 @@ function create_cook_grid(nx, ny)
                Tensors.Vec{2}((0.0,  44.0))]
     grid = generate_grid(Quadrilateral, (nx, ny), corners);
     # facesets for boundary conditions
-    addfaceset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0);
-    addfaceset!(grid, "traction", x -> norm(x[1]) ≈ 48.0);
+    addfacetset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0);
+    addfacetset!(grid, "traction", x -> norm(x[1]) ≈ 48.0);
     return grid
 end;
 
 function create_values(interpolation_u, interpolation_p)
     # quadrature rules
-    qr      = QuadratureRule{2,RefCube}(3)
-    face_qr = QuadratureRule{1,RefCube}(3)
+    qr      = QuadratureRule{RefQuadrilateral}(3)
+    face_qr = FacetQuadratureRule{RefQuadrilateral}(3)
 
     # geometric interpolation
-    interpolation_geom = Lagrange{2,RefCube,1}()
+    interpolation_geom = Lagrange{RefQuadrilateral,1}()
 
     # cell and facevalues for u
-    cellvalues_u = CellVectorValues(qr, interpolation_u, interpolation_geom)
-    facevalues_u = FaceVectorValues(face_qr, interpolation_u, interpolation_geom)
+    cellvalues_u = CellValues(qr, interpolation_u, interpolation_geom)
+    facevalues_u = FacetValues(face_qr, interpolation_u, interpolation_geom)
 
     # cellvalues for p
-    cellvalues_p = CellScalarValues(qr, interpolation_p, interpolation_geom)
+    cellvalues_p = CellValues(qr, interpolation_p, interpolation_geom)
 
     return cellvalues_u, cellvalues_p, facevalues_u
 end;
 
 function create_dofhandler(grid, ipu, ipp)
     dh = DofHandler(grid)
-    push!(dh, :u, 2, ipu) # displacement
-    push!(dh, :p, 1, ipp) # pressure
+    add!(dh, :u, ipu) # displacement
+    add!(dh, :p, ipp) # pressure
     close!(dh)
     return dh
 end;
 
 function create_bc(dh)
     dbc = ConstraintHandler(dh)
-    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "clamped"), (x,t) -> zero(Tensors.Vec{2}), [1,2]))
+    add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "clamped"), (x,t) -> zero(Tensors.Vec{2}), [1,2]))
     close!(dbc)
     t = 0.0
     update!(dbc, t)
@@ -53,28 +53,28 @@ struct LinearElasticity{T}
     K::T
 end
 
-function doassemble(cellvalues_u::CellVectorValues{dim}, cellvalues_p::CellScalarValues{dim},
-                    facevalues_u::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid,
-                    dh::DofHandler, mp::LinearElasticity) where {dim}
+function doassemble(cellvalues_u::CellValues, cellvalues_p::CellValues,
+                    facevalues_u::FacetValues, K::SparseMatrixCSC, grid::Grid{sdim},
+                    dh::DofHandler, mp::LinearElasticity) where {sdim}
 
     f = zeros(ndofs(dh))
     assembler = start_assemble(K, f)
     nu = getnbasefunctions(cellvalues_u)
     np = getnbasefunctions(cellvalues_p)
 
-    fe = PseudoBlockArray(zeros(nu + np), [nu, np]) # local force vector
-    ke = PseudoBlockArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix
+    fe = BlockedArray(zeros(nu + np), [nu, np]) # local force vector
+    ke = BlockedArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix
 
     # traction vector
     t = Tensors.Vec{2}((0.0, 1/16))
     # cache ɛdev outside the element routine to avoid some unnecessary allocations
-    ɛdev = [zero(SymmetricTensor{2, dim}) for i in 1:getnbasefunctions(cellvalues_u)]
+    ɛdev = [zero(SymmetricTensor{2, sdim}) for i in 1:getnbasefunctions(cellvalues_u)]
 
     for cell in CellIterator(dh)
         fill!(ke, 0)
         fill!(fe, 0)
         assemble_up!(ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, grid, mp, ɛdev, t)
-        assemble!(assembler, celldofs(cell), fe, ke)
+        assemble!(assembler, celldofs(cell), ke, fe)
     end
 
     return K, f
@@ -121,8 +121,8 @@ function assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, gr
     # We integrate the Neumann boundary using the facevalues.
     # We loop over all the faces in the cell, then check if the face
     # is in our `"traction"` faceset.
-    @inbounds for face in 1:nfaces(cell)
-        if onboundary(cell, face) && (cellid(cell), face) ∈ getfaceset(grid, "traction")
+    @inbounds for face in 1:nfacets(cell)
+        if (cellid(cell), face) ∈ getfacetset(grid, "traction")
             reinit!(facevalues_u, cell, face)
             for q_point in 1:getnquadpoints(facevalues_u)
                 dΓ = getdetJdV(facevalues_u, q_point)
@@ -154,27 +154,27 @@ function solve(interpolation_u, interpolation_p, mp)
     cellvalues_u, cellvalues_p, facevalues_u = create_values(interpolation_u, interpolation_p)
 
     # assembly and solve
-    K = create_sparsity_pattern(dh);
+    K = allocate_matrix(dh);
     K, f = doassemble(cellvalues_u, cellvalues_p, facevalues_u, K, grid, dh, mp);
     apply!(K, f, dbc)
-    u = Symmetric(K) \ f;
+    u = K \ f;
 
     # export
-    filename = "cook_" * (isa(interpolation_u, Lagrange{2,RefCube,1}) ? "linear" : "quadratic") *
-                         "_linear"
-    vtk_grid(filename, dh) do vtkfile
-        vtk_point_data(vtkfile, dh, u)
-    end
+    # filename = "cook_" * (isa(interpolation_u, Lagrange{RefQuadrilateral,1}) ? "linear" : "quadratic") *
+    #                      "_linear"
+    # vtk_grid(filename, dh) do vtkfile
+    #     vtk_point_data(vtkfile, dh, u)
+    # end
     return u,dh
 end
 
-linear    = Lagrange{2,RefCube,1}()
-quadratic = Lagrange{2,RefCube,2}()
+linear    = Lagrange{RefQuadrilateral,1}()
+quadratic = Lagrange{RefQuadrilateral,2}()
 
 ν = 0.4999999
 Emod = 1.
 Gmod = Emod / 2(1 + ν)
 Kmod = Emod * ν / ((1+ν) * (1-2ν))
 mp = LinearElasticity(Gmod, Kmod)
-u_linear,dh_linear = solve(linear, linear, mp);
-u_quadratic,dh_quadratic = solve(quadratic, linear, mp);
+u_linear,dh_linear = solve(linear^2, linear, mp);
+u_quadratic,dh_quadratic = solve(quadratic^2, linear, mp);
diff --git a/docs/src/ferrite-examples/linear-elasticity.jl b/docs/src/ferrite-examples/linear-elasticity.jl
new file mode 100644
index 0000000..fd8b81c
--- /dev/null
+++ b/docs/src/ferrite-examples/linear-elasticity.jl
@@ -0,0 +1,107 @@
+using Ferrite, SparseArrays
+
+function assemble_linear_elastic_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues, coords::Vector, rhs::Function)
+    n_basefuncs = getnbasefunctions(cellvalues)
+
+    E = 200e3 # Young's modulus [MPa]
+    ν = 0.3 # Poisson's ratio [-]
+
+    G = E/2(1+ν)
+    K = E/3(1-2ν)
+
+    fill!(Ke, 0)
+    fill!(fe, 0)
+
+    for q_point in 1:getnquadpoints(cellvalues)
+        dΩ = getdetJdV(cellvalues, q_point)
+        for i in 1:n_basefuncs
+            Nᵢ = shape_value(cellvalues, q_point, i)
+            x = spatial_coordinate(cellvalues, q_point, coords)
+            fe[i] += rhs(x) ⋅ Nᵢ * dΩ
+            ∇ˢʸᵐNᵢ = shape_symmetric_gradient(cellvalues, q_point, i)
+            for j in 1:n_basefuncs
+                ∇ˢʸᵐNⱼ = shape_symmetric_gradient(cellvalues, q_point, j)
+                σ = 2G * dev(∇ˢʸᵐNⱼ) + K * tr(∇ˢʸᵐNⱼ) * one(∇ˢʸᵐNⱼ)
+                Ke[i, j] += σ ⊡ ∇ˢʸᵐNᵢ * dΩ
+            end
+        end
+    end
+
+    return Ke, fe
+end
+
+function assemble_steady_linear_elastic_global(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler, rhs::Function)
+    n_basefuncs = getnbasefunctions(cellvalues)
+    Ke = zeros(n_basefuncs, n_basefuncs)
+    fe = zeros(n_basefuncs)
+
+    f = zeros(ndofs(dh))
+
+    assembler = start_assemble(K, f)
+
+    for cell in CellIterator(dh)
+        reinit!(cellvalues, cell)
+        assemble_linear_elastic_element!(Ke, fe, cellvalues, getcoordinates(cell), rhs)
+        assemble!(assembler, celldofs(cell), Ke, fe)
+    end
+    return K, f
+end
+
+function manufactured_linear_elastic_problem(element_type, ip, num_elements_per_dim, component)
+    dim = Ferrite.getrefdim(ip)
+    grid = generate_grid(element_type, ntuple(x->num_elements_per_dim, dim));
+    ip_geo = Ferrite.geometric_interpolation(typeof(grid.cells[1]))
+    qr = QuadratureRule{Ferrite.getrefshape(ip)}(2*Ferrite.getorder(ip))
+    cellvalues = CellValues(qr, ip, ip_geo);
+
+    ∂Ω = union(
+        getfacetset(grid, "left"),
+        getfacetset(grid, "right"),
+        getfacetset(grid, "top"),
+        getfacetset(grid, "bottom"),
+    );
+
+    if dim == 3
+        ∂Ω = union(
+            ∂Ω,
+            getfacetset(grid, "front"),
+            getfacetset(grid, "back")
+        )
+    end
+
+    dh = DofHandler(grid)
+    add!(dh, :u, ip)
+    close!(dh);
+
+    K = allocate_matrix(dh)
+
+    ch = ConstraintHandler(dh);
+    dbc = Dirichlet(:u, ∂Ω, (x, t) -> [0.0 for i ∈ 1:dim])
+    add!(ch, dbc);
+
+    close!(ch)
+    update!(ch, 0.0);
+
+    function analytical_rhs(x_eval)
+        E = 200e3 # Young's modulus [MPa]
+        ν = 0.3 # Poisson's ratio [-]
+
+        G = E/2(1+ν)
+        K = E/3(1-2ν)
+
+        u(x)  = Vec{2}(i->i==component ? (π/2)^2 * 2 * prod(cos, x*π/2) : 0.0)
+        ∇u(x) = Tensors.gradient(u, x)
+        ε(x)  = (∇u(x)  + transpose(∇u(x)))/2
+        σ(x)  = 2G * dev(ε(x)) + K * tr(ε(x)) * one(ε(x))
+
+        return Tensors.divergence(σ, x_eval)
+    end
+
+    K, f = assemble_steady_linear_elastic_global(cellvalues, K, dh, analytical_rhs);
+    apply!(K, f, ch)
+    u = K \ f;
+
+    return dh, u
+end
+
+Tensors.tr(t::Tensor{3,dim}) where dim = Vec{dim}(i->sum(t[:,i,:]))
diff --git a/docs/src/ferrite-examples/mixedgrid.jl b/docs/src/ferrite-examples/mixedgrid.jl
index d2f0a6f..0c8ec34 100644
--- a/docs/src/ferrite-examples/mixedgrid.jl
+++ b/docs/src/ferrite-examples/mixedgrid.jl
@@ -1,75 +1,28 @@
-using FerriteGmsh
 using Ferrite
 
-gmsh.initialize()
-gmsh.option.setNumber("General.Terminal", 1)
-gmsh.model.add("demo")
+elements = [Triangle((4, 17, 9)), Triangle((11, 16, 5)), Triangle((5, 16, 12)), Triangle((13, 17, 4)), Triangle((9, 18, 3)), Triangle((3, 18, 11)), Triangle((12, 15, 6)), Triangle((6, 15, 13)), Triangle((17, 18, 9)), Triangle((12, 16, 15)), Triangle((16, 17, 15)), Triangle((11, 18, 16)), Triangle((16, 18, 17)), Triangle((15, 17, 13)), Quadrilateral((1, 7, 14, 10)), Quadrilateral((10, 14, 9, 3)), Quadrilateral((7, 2, 8, 14)), Quadrilateral((14, 8, 4, 9))]
+nodes = Node.(Vec{2, Float64}.([[-0.5, -1.0], [0.5, -1.0], [-0.5, 0.0], [0.5, 0.0], [-0.5, 1.0], [0.5, 1.0], [-1.3751222383007189e-12, -1.0], [0.5, -0.5000000000020595], [1.3751222383007189e-12, 0.0], [-0.5, -0.5000000000020595], [-0.5, 0.4999999999986921], [-1.3751222383007189e-12, 1.0], [0.5, 0.4999999999986921], [-5.664046543035832e-24, -0.5000000000020595], [0.206249999999433, 0.7062499999994473], [-0.12500000000067074, 0.6249999999999999], [0.14375000000014615, 0.35208333333295083], [-0.21874999999926012, 0.28124999999956607]]))
+facesets = Dict{String, Set{FaceIndex}}("bottom" => Set([FaceIndex((15, 1)), FaceIndex((17, 1))]), "top" => Set([FaceIndex((7, 3)), FaceIndex((3, 3))]))
+cellsets = Dict{String, Set{Int64}}("quad" => Set([15, 16, 18, 17]), "triangle" => Set([5, 12, 8, 1, 6, 11, 9, 14, 3, 7, 4, 13, 2, 10]))
 
-lc = 0.2
-gmsh.model.geo.addPoint(-0.5, -1, 0, lc, 1)
-gmsh.model.geo.addPoint(0.5, -1, 0, lc, 2)
-gmsh.model.geo.addPoint(-0.5, 0, 0, lc, 3)
-gmsh.model.geo.addPoint(0.5, 0, 0, lc, 4)
-gmsh.model.geo.addPoint(-0.5, 1, 0, lc, 5)
-gmsh.model.geo.addPoint(0.5, 1, 0, lc, 6)
-
-gmsh.model.geo.addLine(1, 2, 1)
-gmsh.model.geo.addLine(2, 4, 2)
-gmsh.model.geo.addLine(4, 3, 3)
-gmsh.model.geo.addLine(1, 3, 4)
-gmsh.model.geo.addLine(3, 5, 5)
-gmsh.model.geo.addLine(5, 6, 6)
-gmsh.model.geo.addLine(4, 6, 7)
-
-gmsh.model.geo.addCurveLoop([1, 2, 3, -4], 1)
-gmsh.model.geo.addCurveLoop([-3, 7, -6, -5], 2)
-gmsh.model.geo.addPlaneSurface([1], 1)
-gmsh.model.geo.addPlaneSurface([2], 2)
-gmsh.model.geo.mesh.setTransfiniteCurve(1, 3)
-gmsh.model.geo.mesh.setTransfiniteCurve(2, 3)
-gmsh.model.geo.mesh.setTransfiniteCurve(3, 3)
-gmsh.model.geo.mesh.setTransfiniteCurve(4, 3)
-gmsh.model.geo.mesh.setTransfiniteCurve(5, 3)
-gmsh.model.geo.mesh.setTransfiniteCurve(6, 3)
-gmsh.model.geo.mesh.setTransfiniteCurve(7, 3)
-gmsh.model.geo.mesh.setTransfiniteSurface(1)
-gmsh.model.geo.mesh.setRecombine(2, 1)
-
-gmsh.model.addPhysicalGroup(2, [1], 1)
-gmsh.model.setPhysicalName(2, 1, "quad")
-
-gmsh.model.addPhysicalGroup(2, [2], 2)
-gmsh.model.setPhysicalName(2, 2, "triangle")
-
-gmsh.model.addPhysicalGroup(1, [6], 3)
-gmsh.model.setPhysicalName(1, 3, "top")
-
-gmsh.model.addPhysicalGroup(1, [1], 4)
-gmsh.model.setPhysicalName(1, 4, "bottom")
-
-gmsh.model.geo.synchronize()
-gmsh.model.mesh.generate(2)
-
-nodes = tonodes()
-elements, gmsh_eleidx = toelements(2)
-boundarydict = toboundary(1)
-facesets = tofacesets(boundarydict,elements)
-cellsets = tocellsets(2,gmsh_eleidx)
 grid = Grid(elements,nodes,facesets=facesets,cellsets=cellsets)
 
-dh = MixedDofHandler(grid)
-push!(dh,FieldHandler([Field(:p,Lagrange{2,RefTetrahedron,1}(),1),Field(:u,Lagrange{2,RefTetrahedron,1}(),2)], getcellset(grid,"triangle")))
-push!(dh,FieldHandler([Field(:u,Lagrange{2,RefCube,1}(),2)], getcellset(grid,"quad")))
+dh = DofHandler(grid)
+sdh_tri = SubDofHandler(dh, getcellset(grid, "triangle"))
+add!(sdh_tri, :u, Lagrange{RefTriangle,1}()^2)
+add!(sdh_tri, :p, Lagrange{RefTriangle,1}())
+sdh_quad = SubDofHandler(dh, getcellset(grid, "quad"))
+add!(sdh_quad, :u, Lagrange{RefQuadrilateral,1}()^2)
 close!(dh)
 
 u = zeros(ndofs(dh))
 
-for cell in CellIterator(dh,collect(dh.fieldhandlers[1].cellset))
+for cell in CellIterator(dh,collect(dh.subdofhandlers[1].cellset))
     celldofs_ = celldofs(cell)
     u[celldofs_] .= 1
 end
-for cell in CellIterator(dh,collect(dh.fieldhandlers[2].cellset))
+for cell in CellIterator(dh,collect(dh.subdofhandlers[2].cellset))
     celldofs_ = celldofs(cell)
-    dof_range_ = Ferrite.dof_range(dh.fieldhandlers[2],:u)
+    dof_range_ = Ferrite.dof_range(dh.subdofhandlers[2],:u)
     u[celldofs_[dof_range_]] .= 0.5
 end
diff --git a/docs/src/ferrite-examples/plasticity-live.jl b/docs/src/ferrite-examples/plasticity-live.jl
index 99cb794..a24cbdf 100644
--- a/docs/src/ferrite-examples/plasticity-live.jl
+++ b/docs/src/ferrite-examples/plasticity-live.jl
@@ -117,7 +117,7 @@ end;
 function create_dofhandler(grid, interpolation)
     dh = DofHandler(grid)
     dim = 3
-    push!(dh, :u, dim, interpolation) # add a displacement field with 3 components
+    add!(dh, :u, interpolation^3) # add a displacement field with 3 components
     close!(dh)
     return dh
 end
@@ -237,7 +237,7 @@ function solve()
     display(fig)
     Δu = zeros(n_dofs)  # displacement correction
     r = zeros(n_dofs)   # residual
-    K = create_sparsity_pattern(dh); # tangent stiffness matrix
+    K = allocate_matrix(dh); # tangent stiffness matrix
 
     # Create material states. One array for each cell, where each element is an array of material-
     # states - one for each integration point
diff --git a/docs/src/ferrite-examples/plasticity-timeslider.jl b/docs/src/ferrite-examples/plasticity-timeslider.jl
index 73a6cda..7dd17fa 100644
--- a/docs/src/ferrite-examples/plasticity-timeslider.jl
+++ b/docs/src/ferrite-examples/plasticity-timeslider.jl
@@ -235,7 +235,7 @@ function solve()
     u_history = Vector{Vector{Float64}}()
     Δu = zeros(n_dofs)  # displacement correction
     r = zeros(n_dofs)   # residual
-    K = create_sparsity_pattern(dh); # tangent stiffness matrix
+    K = allocate_matrix(dh); # tangent stiffness matrix
 
     # Create material states. One array for each cell, where each element is an array of material-
     # states - one for each integration point
diff --git a/docs/src/ferrite-examples/plasticity.jl b/docs/src/ferrite-examples/plasticity.jl
index e5cf900..62102f5 100644
--- a/docs/src/ferrite-examples/plasticity.jl
+++ b/docs/src/ferrite-examples/plasticity.jl
@@ -101,23 +101,22 @@ end
 
 function create_values(interpolation)
     # setup quadrature rules
-    qr      = QuadratureRule{3,RefTetrahedron}(2)
-    face_qr = QuadratureRule{2,RefTetrahedron}(3)
+    qr      = QuadratureRule{RefTetrahedron}(2)
+    face_qr = FacetQuadratureRule{RefTetrahedron}(3)
 
     # create geometric interpolation (use the same as for u)
-    interpolation_geom = Lagrange{3,RefTetrahedron,1}()
+    interpolation_geom = Lagrange{RefTetrahedron,1}()^3
 
     # cell and facevalues for u
-    cellvalues_u = CellVectorValues(qr, interpolation, interpolation_geom)
-    facevalues_u = FaceVectorValues(face_qr, interpolation, interpolation_geom)
+    cellvalues_u = CellValues(qr, interpolation, interpolation_geom)
+    facevalues_u = FacetValues(face_qr, interpolation, interpolation_geom)
 
     return cellvalues_u, facevalues_u
 end;
 
 function create_dofhandler(grid, interpolation)
     dh = DofHandler(grid)
-    dim = 3
-    push!(dh, :u, dim, interpolation) # add a displacement field with 3 components
+    add!(dh, :u, interpolation) # add a displacement field with 3 components
     close!(dh)
     return dh
 end
@@ -126,14 +125,14 @@ function create_bc(dh, grid)
     dbcs = ConstraintHandler(dh)
     # Clamped on the left side
     dofs = [1, 2, 3]
-    dbc = Dirichlet(:u, getfaceset(grid, "left"), (x,t) -> [0.0, 0.0, 0.0], dofs)
+    dbc = Dirichlet(:u, getfacetset(grid, "left"), (x,t) -> [0.0, 0.0, 0.0], dofs)
     add!(dbcs, dbc)
     close!(dbcs)
     return dbcs
 end;
 
-function doassemble(cellvalues::CellVectorValues{dim},
-                    facevalues::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid,
+function doassemble(cellvalues::CellValues{dim},
+                    facevalues::FacetValues{dim}, K::SparseMatrixCSC, grid::Grid,
                     dh::DofHandler, material::J2Plasticity, u, states, t) where {dim}
     r = zeros(ndofs(dh))
     assembler = start_assemble(K, r)
@@ -148,7 +147,7 @@ function doassemble(cellvalues::CellVectorValues{dim},
         ue = u[eldofs]
         assemble_cell!(ke, re, cell, cellvalues, facevalues, grid, material,
                        ue, state, t)
-        assemble!(assembler, eldofs, re, ke)
+        assemble!(assembler, eldofs, ke, re)
     end
     return K, r
 end
@@ -178,8 +177,8 @@ function assemble_cell!(Ke, re, cell, cellvalues, facevalues, grid, material,
     symmetrize_lower!(Ke)
 
     # Add traction as a negative contribution to the element residual `re`:
-    for face in 1:nfaces(cell)
-        if onboundary(cell, face) && (cellid(cell), face) ∈ getfaceset(grid, "right")
+    for face in 1:nfacets(cell)
+        if (cellid(cell), face) ∈ getfacetset(grid, "right")
             reinit!(facevalues, cell, face)
             for q_point in 1:getnquadpoints(facevalues)
                 dΓ = getdetJdV(facevalues, q_point)
@@ -219,15 +218,15 @@ function solve(liveplotting=false)
     # Create geometry, dofs and boundary conditions
     n = 2
     nels = (10n, n, 2n) # number of elements in each spatial direction
-    P1 = Vec((0.0, 0.0, 0.0))  # start point for geometry
-    P2 = Vec((L, w, h))        # end point for geometry
+    P1 = Tensors.Vec((0.0, 0.0, 0.0))  # start point for geometry
+    P2 = Tensors.Vec((L, w, h))        # end point for geometry
     grid = generate_grid(Tetrahedron, nels, P1, P2)
-    interpolation = Lagrange{3, RefTetrahedron, 1}() # Linear tet with 3 unknowns/node
+    interpolation = Lagrange{RefTetrahedron, 1}() # Linear tet with 3 unknowns/node
 
-    dh = create_dofhandler(grid, interpolation) # JuaFEM helper function
+    dh = create_dofhandler(grid, interpolation^3) # helper function
     dbcs = create_bc(dh, grid) # create Dirichlet boundary-conditions
 
-    cellvalues, facevalues = create_values(interpolation)
+    cellvalues, facevalues = create_values(interpolation^3)
 
     # Pre-allocate solution vectors, etc.
     n_dofs = ndofs(dh)  # total number of dofs
@@ -240,7 +239,7 @@ function solve(liveplotting=false)
     end
     Δu = zeros(n_dofs)  # displacement correction
     r = zeros(n_dofs)   # residual
-    K = create_sparsity_pattern(dh); # tangent stiffness matrix
+    K = allocate_matrix(dh); # tangent stiffness matrix
 
     # Create material states. One array for each cell, where each element is an array of material-
     # states - one for each integration point
@@ -255,7 +254,7 @@ function solve(liveplotting=false)
 
     for timestep in 1:n_timesteps
         t = timestep # actual time (used for evaluating d-bndc)
-        traction = Vec((0.0, 0.0, traction_magnitude[timestep]))
+        traction = Tensors.Vec((0.0, 0.0, traction_magnitude[timestep]))
         newton_itr = -1
         update!(dbcs, t) # evaluates the D-bndc at time t
         apply!(u, dbcs)  # set the prescribed values in the solution vector
diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md
index 98f1812..a67fcf9 100644
--- a/docs/src/tutorial.md
+++ b/docs/src/tutorial.md
@@ -15,8 +15,10 @@ and solution vector because we need to pass those objects to `MakiePlotter`.
 ### Mesh utilities
 
 ```@example 1
-import JSServe # hide
-JSServe.Page() # hide
+import WGLMakie, Bonito # hide
+Bonito.Page() # hide
+WGLMakie.activate!() # hide
+WGLMakie.Makie.inline!(true) # hide
 ```
 
 You can start by plotting your mesh
@@ -25,10 +27,10 @@ You can start by plotting your mesh
 import FerriteViz
 using Ferrite
 import WGLMakie #activating the backend, switch to GLMakie or CairoMakie (for 2D) locally
-WGLMakie.set_theme!(resolution=(800, 400)) # hide
+WGLMakie.set_theme!(size=(800, 400)) # hide
 
 grid = generate_grid(Hexahedron,(3,3,3))
-FerriteViz.wireframe(grid,markersize=10,strokewidth=2)
+FerriteViz.wireframe(grid,markersize=10,linewidth=2)
 ```
 
 FerriteViz.jl also supports showing labels for `Ferrite.AbstractGrid` entities, such as node- and celllabels, as well as plotting cellsets.
@@ -38,7 +40,7 @@ grid = generate_grid(Quadrilateral,(3,3))
 addcellset!(grid,"s1",Set((1,4,7)))
 addcellset!(grid,"s2",Set((2,5,8)))
 addcellset!(grid,"s3",Set((3,6,9)))
-FerriteViz.wireframe(grid,markersize=10,strokewidth=1,nodelabels=true,celllabels=true,cellsets=true)
+FerriteViz.wireframe(grid,markersize=10,linewidth=1,nodelabels=true,celllabels=true,cellsets=true)
 ```
 
 ### Solution field of a boundary value problem
@@ -80,7 +82,7 @@ u, dh, uhistory, σ, κ = solve()
 plotter = FerriteViz.MakiePlotter(dh,u)
 
 FerriteViz.solutionplot(plotter,colormap=:thermal,deformation_field=:u)
-FerriteViz.wireframe!(plotter,deformation_field=:u,markersize=10,strokewidth=1)
+FerriteViz.wireframe!(plotter,deformation_field=:u,markersize=10,linewidth=1)
 WGLMakie.current_figure()
 ```
 
@@ -90,7 +92,7 @@ FerriteViz.jl also supports to plot cell data, such as the **averaged** von-Mise
 ```@example 1
 u, dh, uhistory, σ, κ = solve()
 FerriteViz.cellplot(plotter,σ,colormap=:thermal,deformation_field=:u,deformation_scale=2.0)
-FerriteViz.wireframe!(plotter,deformation_field=:u,markersize=10,strokewidth=1,deformation_scale=2.0)
+FerriteViz.wireframe!(plotter,deformation_field=:u,markersize=10,linewidth=1,deformation_scale=2.0)
 WGLMakie.current_figure()
 ```
 For a more granular investigation of the stress field consult the advanced tutorial.
diff --git a/src/FerriteViz.jl b/src/FerriteViz.jl
index 3a35828..846ab17 100644
--- a/src/FerriteViz.jl
+++ b/src/FerriteViz.jl
@@ -1,8 +1,10 @@
 module FerriteViz
 
 using Makie
+using StaticArrays
 using Tensors
-import Ferrite
+using Ferrite
+import Ferrite: get_grid, getrefshape
 import GeometryBasics
 import ShaderAbstractions
 import LinearAlgebra
diff --git a/src/lor_tools.jl b/src/lor_tools.jl
index 7a20503..a7c30fa 100644
--- a/src/lor_tools.jl
+++ b/src/lor_tools.jl
@@ -1,17 +1,17 @@
 # These functions generate the corresponding first order cells of an interpolation.
 # Triangle
-for_nodes(::Union{Ferrite.Lagrange{2,Ferrite.RefTetrahedron,1},Ferrite.DiscontinuousLagrange{2,Ferrite.RefTetrahedron,1},Ferrite.Triangle}) = (
+for_nodes(::Union{Lagrange{RefTriangle,1},DiscontinuousLagrange{RefTriangle,1},Triangle}) = (
     (3,1,2),
 )
 # Quadratic Triangle
-for_nodes(::Union{Ferrite.Lagrange{2,Ferrite.RefTetrahedron,2},Ferrite.DiscontinuousLagrange{2,Ferrite.RefTetrahedron,2},Ferrite.QuadraticTriangle}) = (
+for_nodes(::Union{Lagrange{RefTriangle,2},DiscontinuousLagrange{RefTriangle,2},QuadraticTriangle}) = (
     (6,1,4),
     (5,6,4),
     (3,6,5),
     (5,4,2),
 )
 # Cubic Triangle
-for_nodes(::Union{Ferrite.Lagrange{2,Ferrite.RefTetrahedron,3},Ferrite.DiscontinuousLagrange{2,Ferrite.RefTetrahedron,3},Ferrite.Cell{2,10,3}}) = (
+for_nodes(::Union{Lagrange{RefTriangle,3},DiscontinuousLagrange{RefTriangle,3}}) = (
     (3,8,7),
     (7,8,10),
     (8,9,10),
@@ -23,7 +23,7 @@ for_nodes(::Union{Ferrite.Lagrange{2,Ferrite.RefTetrahedron,3},Ferrite.Discontin
     (10,4,5),
 )
 # Biquadratic Triangle
-for_nodes(::Union{Ferrite.Lagrange{2,Ferrite.RefTetrahedron,4},Ferrite.DiscontinuousLagrange{2,Ferrite.RefTetrahedron,4},Ferrite.Cell{2,15,3}}) = (
+for_nodes(::Union{Lagrange{RefTriangle,4},DiscontinuousLagrange{RefTriangle,4}}) = (
     (3,10,9),
     (13,9,10),
     (10,11,13),
@@ -42,7 +42,7 @@ for_nodes(::Union{Ferrite.Lagrange{2,Ferrite.RefTetrahedron,4},Ferrite.Discontin
     (7,6,2),
 )
 # Quintic Triangle
-for_nodes(::Union{Ferrite.Lagrange{2,Ferrite.RefTetrahedron,5},Ferrite.DiscontinuousLagrange{2,Ferrite.RefTetrahedron,5},Ferrite.Cell{2,20,3}}) = (
+for_nodes(::Union{Lagrange{RefTriangle,5},DiscontinuousLagrange{RefTriangle,5}}) = (
     (3,12,11),
     (16,11,12),
     (12,13,16),
@@ -70,11 +70,11 @@ for_nodes(::Union{Ferrite.Lagrange{2,Ferrite.RefTetrahedron,5},Ferrite.Discontin
     (8,7,2),
 )
 # Tetrahedron
-for_nodes(::Union{Ferrite.Lagrange{3,Ferrite.RefTetrahedron,1},Ferrite.DiscontinuousLagrange{3,Ferrite.RefTetrahedron,1},Ferrite.Tetrahedron}) = (
+for_nodes(::Union{Lagrange{RefTetrahedron,1},DiscontinuousLagrange{RefTetrahedron,1},Tetrahedron}) = (
     (1,2,3,4),
 )
 # Quadratic Tetrahedron
-for_nodes(::Union{Ferrite.Lagrange{3,Ferrite.RefTetrahedron,2},Ferrite.DiscontinuousLagrange{3,Ferrite.RefTetrahedron,2},Ferrite.QuadraticTetrahedron}) = (
+for_nodes(::Union{Lagrange{RefTetrahedron,2},DiscontinuousLagrange{RefTetrahedron,2},QuadraticTetrahedron}) = (
     (5,2,6,9),
     (7,6,3,10),
     (8,9,10,4),
@@ -85,22 +85,22 @@ for_nodes(::Union{Ferrite.Lagrange{3,Ferrite.RefTetrahedron,2},Ferrite.Discontin
     (9,10,8,6),
 )
 # Quadrilateral
-for_nodes(::Union{Ferrite.Lagrange{2,Ferrite.RefCube,1},Ferrite.DiscontinuousLagrange{2,Ferrite.RefCube,1},Ferrite.Quadrilateral}) = (
+for_nodes(::Union{Lagrange{RefQuadrilateral,1},DiscontinuousLagrange{RefQuadrilateral,1},Quadrilateral}) = (
     (1,2,3,4),
 )
 # Quadratic Quadrilateral
-for_nodes(::Union{Ferrite.Lagrange{2,Ferrite.RefCube,2},Ferrite.DiscontinuousLagrange{2,Ferrite.RefCube,2},Ferrite.QuadraticQuadrilateral}) = (
+for_nodes(::Union{Lagrange{RefQuadrilateral,2},DiscontinuousLagrange{RefQuadrilateral,2},QuadraticQuadrilateral}) = (
     (1,5,9,8),
     (5,2,6,9),
     (9,6,3,7),
     (8,9,7,4),
 )
 # Hexahedron
-for_nodes(::Union{Ferrite.Lagrange{3,Ferrite.RefCube,1},Ferrite.DiscontinuousLagrange{3,Ferrite.RefCube,1},Ferrite.Hexahedron}) = (
+for_nodes(::Union{Lagrange{RefHexahedron,1},DiscontinuousLagrange{RefHexahedron,1},Hexahedron}) = (
     (1,2,3,4,5,6,7,8),
 )
 # Quadratic Hexahedron
-for_nodes(::Union{Ferrite.Lagrange{3,Ferrite.RefCube,2},Ferrite.DiscontinuousLagrange{3,Ferrite.RefCube,2},Ferrite.Cell{3,27,6}}) = (
+for_nodes(::Union{Lagrange{RefHexahedron,2},DiscontinuousLagrange{RefHexahedron,2},QuadraticHexahedron}) = (
     (1,9,21,12,17,22,27,25),
     (17,22,27,25,5,13,26,16),
     (9,2,10,21,22,18,23,27),
@@ -114,21 +114,12 @@ for_nodes(::Union{Ferrite.Lagrange{3,Ferrite.RefCube,2},Ferrite.DiscontinuousLag
 """
 Get the interpolation of the first order refinement. 
 """
-for_interpolation(ip::Ferrite.Lagrange{dim,shape,order}) where {dim,shape,order} = Ferrite.Lagrange{dim,shape,1}()
+for_interpolation(ip::Lagrange{shape,order}) where {shape,order} = Lagrange{shape,1}()
 
-for_base_geometry_type(ip::Ferrite.Lagrange{2,Ferrite.RefCube,order}) where {order} = Ferrite.Quadrilateral
-for_base_geometry_type(ip::Ferrite.Lagrange{3,Ferrite.RefCube,order}) where {order} = Ferrite.Hexahedron
-for_base_geometry_type(ip::Ferrite.Lagrange{2,Ferrite.RefTetrahedron,order}) where {order} = Ferrite.Triangle
-for_base_geometry_type(ip::Ferrite.Lagrange{3,Ferrite.RefTetrahedron,order}) where {order} = Ferrite.Tetrahedron
-
-# TODO move into ferrite core
-function Ferrite.field_offset(dh::Ferrite.DofHandler, field_name::Int)
-    offset = 0
-    for i in 1:field_name-1
-        offset += Ferrite.getnbasefunctions(dh.field_interpolations[i])::Int * dh.field_dims[i]
-    end
-    return offset
-end
+for_base_geometry_type(ip::Lagrange{RefQuadrilateral,order}) where {order} = Quadrilateral
+for_base_geometry_type(ip::Lagrange{RefHexahedron,order}) where {order} = Hexahedron
+for_base_geometry_type(ip::Lagrange{RefTriangle,order}) where {order} = Triangle
+for_base_geometry_type(ip::Lagrange{RefTetrahedron,order}) where {order} = Tetrahedron
 
 # TODO move into ferrite core
 function Ferrite.dof_range(dh::Ferrite.DofHandler, field_idx::Int)
@@ -145,20 +136,23 @@ Create a first order discretization w.r.t. a field and transfer
 the solution.
 """
 function for_discretization(dh, u)
+    @assert length(dh.subdofhandlers) == 1 "Subdomains not supported yet"
+    # TODO Dofs for fields are not continuous. Think harder.
+    @assert length(Ferrite.getfieldnames(dh)) == 1 "Multiple fields not supported yet"
+    sdh = dh.subdofhandlers[1]
     field_idx=1
+    grid = dh.grid
 
     # Some helpers
-    ip = Ferrite.getfieldinterpolation(dh, field_idx)
-    field_dim = Ferrite.getfielddim(dh, field_idx)
-    spatial_dim = Ferrite.getdim(dh.grid)
+    ip = Ferrite.getfieldinterpolation(sdh, field_idx)
+    field_dim = Ferrite.n_components(sdh, field_idx)
+    spatial_dim = Ferrite.getspatialdim(grid)
 
-    # TODO Dofs for fields are not continuous. Think harder.
-    @assert Ferrite.nfields(dh) == 1 "Multiple fields not supported yet"
     # # Get dof range, the hard way
     # dof_min = dh.ndofs.x
     # dof_max = 0
     ncells = 0
-    for cell ∈ Ferrite.CellIterator(dh)
+    for cell ∈ Ferrite.CellIterator(sdh)
         # celldofs = Ferrite.celldofs(cell)
         # dof_max = max(dof_max, maximum(celldofs))
         # dof_min = min(dof_min, minimum(celldofs))
@@ -166,21 +160,21 @@ function for_discretization(dh, u)
     end
 
     # Preallocate
-    nodes = Vector{typeof(dh.grid.nodes[1])}(undef, Ferrite.ndofs(dh)) #(dof_max-dof_min+1)÷field_dim)
-    cells = Vector{Ferrite.getcelltype(dh.grid)}(undef, ncells)
+    nodes = Vector{typeof(grid.nodes[1])}(undef, Ferrite.ndofs(dh)) #(dof_max-dof_min+1)÷field_dim)
+    cells = Vector{Ferrite.getcelltype(grid)}(undef, ncells)
 
     ref_coords = Ferrite.reference_coordinates(ip)
     # Starting here we assume a single type of cell being present
     # TODO improve this.
-    ip_geo = Ferrite.default_interpolation(typeof(dh.grid.cells[1]))
+    ip_geo = Ferrite.geometric_interpolation(typeof(grid.cells[1]))
     nodes_per_cell = length(ref_coords)
-    qr = Ferrite.QuadratureRule{spatial_dim, Ferrite.getrefshape(ip)}(zeros(nodes_per_cell), ref_coords)
-    cv = Ferrite.CellScalarValues(qr, ip, ip_geo)
+    qr = Ferrite.QuadratureRule{Ferrite.getrefshape(ip)}(zeros(nodes_per_cell), ref_coords)
+    cv = Ferrite.CellValues(qr, ip, ip_geo)
     cellidx = 1
-    for cell ∈ Ferrite.CellIterator(dh)
+    for cell ∈ Ferrite.CellIterator(sdh)
         Ferrite.reinit!(cv, cell)
         coords = Ferrite.getcoordinates(cell)
-        dofs_f = Ferrite.celldofs(cell)[Ferrite.dof_range(dh, field_idx)]
+        dofs_f = Ferrite.celldofs(cell)[Ferrite.dof_range(sdh, field_idx)]
 
         # Extract coordinates
         for q ∈ 1:nodes_per_cell
@@ -195,17 +189,24 @@ function for_discretization(dh, u)
     end
 
     # Generate a new dof handler.
-    grid_new = Ferrite.Grid(cells, nodes)
-    dh_new = Ferrite.DofHandler(grid_new)
-    Ferrite.push!(dh_new, getfieldname(dh, field_idx), Ferrite.getfielddim(dh, field_idx), for_interpolation(ip))
-    Ferrite.close!(dh_new);
+    grid_new = Grid(cells, nodes)
+    dh_new = DofHandler(grid_new)
+    vdim = Ferrite.n_components(sdh, field_idx)
+    lip = if vdim > 1
+        for_interpolation(ip)^vdim
+    else
+        for_interpolation(ip)
+    end
+    add!(dh_new, getfieldname(sdh, field_idx), lip)
+    close!(dh_new);
 
     # Transfer solution the dumb way.
     # TODO this can be optimized.
     u_new = zeros(Ferrite.ndofs(dh_new))
+    dh_dof_range = Ferrite.celldofs(dh_new, 1)
     for cell_idx ∈ 1:length(dh_new.grid.cells)
-        dh_dof_range = dh_new.cell_dofs_offset[cell_idx]:(dh_new.cell_dofs_offset[cell_idx+1]-1)
-        dofs = dh_new.cell_dofs[dh_dof_range][Ferrite.dof_range(dh_new, field_idx)]
+        Ferrite.celldofs!(dh_dof_range, dh_new, cell_idx)
+        dofs = dh_dof_range[Ferrite.dof_range(dh_new.subdofhandlers[1], 1)]
         u_new[dofs] .= u[[dh_new.grid.cells[cell_idx].nodes...]]
     end
 
diff --git a/src/makieplotting.jl b/src/makieplotting.jl
index dd7daeb..dda53e2 100644
--- a/src/makieplotting.jl
+++ b/src/makieplotting.jl
@@ -22,21 +22,17 @@ keyword arguments are:
 - `colormap::Symbol=:cividis`
 - `colorrange::NTuple{2,<:Number}`: Specify (min, max) of the colorscale. If not given, min and max are calculated automatically from the data. 
 - `deformation_scale=1.0`
-- `shading=false`
-- `scale_plot=false`
-- `transparent=false`
+- `shading=Makie.NoShading`
 - `nan_color::Union{Symbol, <:Colorant}=:red`
 """
 @recipe(SolutionPlot) do scene
     Attributes(
-    scale_plot=false,
-    shading=false,
+    shading=Makie.NoShading,
     field=:default,
     deformation_field=:default,
     process=postprocess,
     colormap=:cividis,
     colorrange=Makie.automatic,
-    transparent=false,
     deformation_scale = 1.0,
     nan_color=:red,
     )
@@ -53,11 +49,13 @@ function Makie.plot!(SP::SolutionPlot{<:Tuple{<:MakiePlotter}})
         end
     end
     u_matrix = @lift begin
+        grid = Ferrite.get_grid(plotter.dh)
+        sdim = Ferrite.getspatialdim(grid)
         if $(SP[:deformation_field])===:default
-             Ferrite.getdim(plotter.dh.grid) > 2 ? Point3f[Point3f(0,0,0)] : Point2f[Point2f(0,0)]
+            sdim > 2 ? Point3f[Point3f(0,0,0)] : Point2f[Point2f(0,0)]
         else
             #TODO remove convert
-            convert(Vector{Point{Ferrite.getdim(plotter.dh.grid),Float32}},Makie.to_vertices(transfer_solution(plotter,$(plotter.u); field_name=$(SP[:deformation_field]), process=identity)))
+            convert(Vector{Point{sdim,Float32}},Makie.to_vertices(transfer_solution(plotter,$(plotter.u); field_name=$(SP[:deformation_field]), process=identity)))
         end
     end
     @lift begin
@@ -67,7 +65,7 @@ function Makie.plot!(SP::SolutionPlot{<:Tuple{<:MakiePlotter}})
             plotter.physical_coords_mesh[1:end] = plotter.physical_coords .+ ($(SP[:deformation_scale]) .* $(u_matrix))
         end
     end
-    return Makie.mesh!(SP, plotter.mesh, color=solution, shading=SP[:shading], scale_plot=SP[:scale_plot], colormap=SP[:colormap],colorrange=SP[:colorrange] , transparent=SP[:transparent], nan_color=SP[:nan_color])
+    return Makie.mesh!(SP, plotter.mesh, color=solution, shading=SP[:shading], colormap=SP[:colormap],colorrange=SP[:colorrange], nan_color=SP[:nan_color])
 end
 
 """
@@ -83,20 +81,16 @@ keyword arguments are:
 - `colormap::Symbol=:cividis`
 - `colorrange::NTuple{2,<:Number}`: Specify (min, max) of the colorscale. If not given, min and max are calculated automatically from the data. 
 - `deformation_scale=1.0`
-- `shading=false`
-- `scale_plot=false`
-- `transparent=false`
+- `shading=Makie.NoShading`
 - `nan_color::Union{Symbol, <:Colorant}=:red`
 """
 @recipe(CellPlot) do scene
     Attributes(
-    scale_plot=false,
-    shading=false,
+    shading=Makie.NoShading,
     deformation_field=:default,
     process=identity,
     colormap=:cividis,
     colorrange=Makie.automatic,
-    transparent=false,
     deformation_scale = 1.0,
     nan_color=:red,
     )
@@ -109,7 +103,9 @@ function Makie.plot!(CP::CellPlot{<:Tuple{<:MakiePlotter{dim},Vector}}) where di
         if $(CP[:deformation_field])===:default
             Point3f[Point3f(0,0,0)]
         else
-            convert(Vector{Point{Ferrite.getdim(plotter.dh.grid),Float32}},Makie.to_vertices(transfer_solution(plotter,$(plotter.u); field_name=$(CP[:deformation_field]), process=identity)))
+            grid = Ferrite.get_grid(plotter.dh)
+            sdim = Ferrite.getspatialdim(grid)
+            convert(Vector{Point{sdim,Float32}},Makie.to_vertices(transfer_solution(plotter,$(plotter.u); field_name=$(CP[:deformation_field]), process=identity)))
         end
     end
     coords = @lift begin
@@ -120,7 +116,7 @@ function Makie.plot!(CP::CellPlot{<:Tuple{<:MakiePlotter{dim},Vector}}) where di
         end
     end
     solution =  @lift(reshape(transfer_scalar_celldata(plotter, qp_values; process=$(CP[:process])), num_vertices(plotter)))
-    return Makie.mesh!(CP, plotter.mesh, color=solution, shading=CP[:shading], scale_plot=CP[:scale_plot], colormap=CP[:colormap], transparent=CP[:transparent], colorrange=CP[:colorrange], nan_color=CP[:nan_color])
+    return Makie.mesh!(CP, plotter.mesh, color=solution, shading=CP[:shading], colormap=CP[:colormap], colorrange=CP[:colorrange], nan_color=CP[:nan_color])
 end
 
 """
@@ -133,7 +129,7 @@ end
 Plots the finite element mesh, optionally labels it and transforms it if a suitable `deformation_field` is given.
 
 - `plotnodes::Bool=true` plots the nodes as circles/spheres
-- `strokewidth::Int=2` how thick faces/edges are drawn
+- `linewidth::Int=2` how thick faces/edges are drawn
 - `color::Symbol=theme(scene,:linecolor)` color of the faces/edges and nodes
 - `markersize::Int=30` size of the nodes
 - `deformation_field::Symbol=:default` field that transforms the mesh by the given deformation, defaults to no deformation
@@ -143,19 +139,19 @@ Plots the finite element mesh, optionally labels it and transforms it if a suita
 - `nodelabelcolor=:darkblue`
 - `celllabels=false` global cell id labels
 - `celllabelcolor=:darkred`
-- `textsize::Int=15` size of the label's text
+- `fontsize::Int=15` size of the label's text
 - `visible=true`
 """
 @recipe(Wireframe) do scene
     Attributes(
     plotnodes=true,
     color=theme(scene, :linecolor),
-    strokewidth=theme(scene, :linewidth),
+    linewidth=theme(scene, :linewidth),
     markersize=theme(scene, :markersize),
     deformation_field=:default,
     visible=true,
     deformation_scale=1,
-    textsize=15,
+    fontsize=15,
     offset=(0.0,0.0),
     nodelabels=false,
     nodelabelcolor=:darkblue,
@@ -179,7 +175,10 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:MakiePlotter{dim}}}) where dim
         if $(WF[:deformation_field])===:default
             pointtype[zero(pointtype)]
         else
-            convert(Vector{Point{Ferrite.getdim(plotter.dh.grid),Float32}},Makie.to_vertices(dof_to_node(plotter.dh, $(WF[1][].u); field_name=$(WF[:deformation_field]))))
+            grid = Ferrite.get_grid(plotter.dh)
+            sdim = Ferrite.getspatialdim(grid)
+            # Float32 to be performant on a larger range of GPUs
+            convert(Vector{Point{sdim,Float32}},Makie.to_vertices(dof_to_node(plotter.dh, $(WF[1][].u); field_name=$(WF[:deformation_field]))))
         end
     end
     gridnodes = @lift begin
@@ -191,21 +190,25 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:MakiePlotter{dim}}}) where dim
     end
     lines = @lift begin
         dim > 2 ? (lines = Point3f[]) : (lines = Point2f[])
-        for cell in Ferrite.getcells(plotter.dh.grid)
-            boundaryentities = dim < 3 ? Ferrite.faces(cell) : Ferrite.edges(cell)
+        grid = Ferrite.get_grid(plotter.dh)
+        for cell in Ferrite.getcells(grid)
+            boundaryentities = Ferrite.edges(cell)
             append!(lines, [$gridnodes[e] for boundary in boundaryentities for e in boundary])
         end
         lines
     end
     nodes = @lift($(WF[:plotnodes]) ? $(gridnodes) : pointtype[zero(pointtype)])
     #plot cellsets
-    cellsets = plotter.dh.grid.cellsets
+    grid = Ferrite.get_grid(plotter.dh)
+    cellsets = grid.cellsets
     cellset_to_value = Dict{String,Int}()
     for (cellsetidx,(cellsetname,cellset)) in enumerate(cellsets)
         cellset_to_value[cellsetname] = cellsetidx
     end
-    cellset_u = zeros(Ferrite.getncells(plotter.dh.grid))
-    for (cellidx,cell) in enumerate(Ferrite.getcells(plotter.dh.grid))
+    cellset_u = zeros(Ferrite.getncells(grid))
+    allcells = Ferrite.getcells(grid)
+    ncells = Ferrite.getncells(grid)
+    for (cellidx,cell) in enumerate(allcells)
         for (cellsetname,cellsetvalue) in cellset_to_value
             if cellidx in cellsets[cellsetname]
                 cellset_u[cellidx] = cellsetvalue
@@ -228,7 +231,7 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:MakiePlotter{dim}}}) where dim
     end
     colorrange = isempty(cellset_to_value) ? (0,1) : (0,maximum(values(cellset_to_value)))
     cellset_u =  reshape(transfer_scalar_celldata(plotter, cellset_u; process=identity), num_vertices(plotter))
-    Makie.mesh!(WF, plotter.mesh, color=cellset_u, shading=false, scale_plot=false, colormap=:darktest, visible=WF[:cellsets])
+    Makie.mesh!(WF, plotter.mesh, color=cellset_u, shading=Makie.NoShading, colormap=:darktest, visible=WF[:cellsets])
     #plot the nodes
     shouldplot = @lift ($(WF[:visible]) && $(WF[:plotnodes]))
     Makie.scatter!(WF,gridnodes,markersize=WF[:markersize], color=WF[:color], visible=shouldplot)
@@ -236,22 +239,26 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:MakiePlotter{dim}}}) where dim
     nodelabels = @lift $(WF[:nodelabels]) ? ["$i" for i in 1:size($gridnodes,1)] : [""]
     nodepositions = @lift $(WF[:nodelabels]) ? $gridnodes : (dim < 3 ? Point2f[Point2f((0,0))] : Point3f[Point3f((0,0,0))])
     #set up celllabels
-    celllabels = @lift $(WF[:celllabels]) ? ["$i" for i in 1:Ferrite.getncells(plotter.dh.grid)] : [""]
-    cellpositions = @lift $(WF[:celllabels]) ? [midpoint(cell,$gridnodes) for cell in Ferrite.getcells(plotter.dh.grid)] : (dim < 3 ? [Point2f((0,0))] : [Point3f((0,0,0))])
-    Makie.text!(WF,nodepositions, text=nodelabels, textsize=WF[:textsize], offset=WF[:offset],color=WF[:nodelabelcolor])
-    Makie.text!(WF,celllabels, position=cellpositions, textsize=WF[:textsize], color=WF[:celllabelcolor], align=(:center,:center))
+    celllabels = @lift $(WF[:celllabels]) ? ["$i" for i in 1:ncells] : [""]
+    cellpositions = @lift $(WF[:celllabels]) ? [midpoint(cell,$gridnodes) for cell in allcells] : (dim < 3 ? [Point2f((0,0))] : [Point3f((0,0,0))])
+    Makie.text!(WF,nodepositions, text=nodelabels, fontsize=WF[:fontsize], offset=WF[:offset],color=WF[:nodelabelcolor])
+    Makie.text!(WF,celllabels, position=cellpositions, fontsize=WF[:fontsize], color=WF[:celllabelcolor], align=(:center,:center))
     #plot edges (3D) /faces (2D) of the mesh
-    Makie.linesegments!(WF,lines,color=WF[:color], linewidth=WF[:strokewidth], visible=WF[:visible], depth_shift=WF[:depth_shift])
+    @show lines, "-----------------------------"
+    Makie.linesegments!(WF,lines,color=WF[:color], linewidth=WF[:linewidth], visible=WF[:visible], depth_shift=WF[:depth_shift])
 end
 
 
 function Makie.plot!(WF::Wireframe{<:Tuple{<:Ferrite.AbstractGrid{dim}}}) where dim
-    grid = WF[1][]
-    coords = [Ferrite.getcoordinates(node)[i] for node in Ferrite.getnodes(grid), i in 1:dim] 
+    @info WF
+    grid   = WF[1][]
+    coords = [Ferrite.get_node_coordinate(node)[i] for node in Ferrite.getnodes(grid), i in 1:dim] 
     coords = Makie.to_vertices(coords)
     dim > 2 ? (lines = Point3f[]) : (lines = Point2f[])
-    for cell in Ferrite.getcells(grid)
-        boundaryentities = dim < 3 ? Ferrite.faces(cell) : Ferrite.edges(cell)
+    allcells = Ferrite.getcells(grid)
+    ncells   = Ferrite.getncells(grid)
+    for cell in allcells
+        boundaryentities = Ferrite.edges(cell)
         append!(lines, [coords[e] for boundary in boundaryentities for e in boundary])
     end
     nodes = @lift($(WF[:plotnodes]) ? coords : Point3f[Point3f(0,0,0)])
@@ -259,21 +266,18 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:Ferrite.AbstractGrid{dim}}}) where
     Makie.scatter!(WF,nodes,markersize=WF[:markersize], color=WF[:color], visible=shouldplot)
     nodelabels = @lift $(WF[:nodelabels]) ? ["$i" for i in 1:size(coords,1)] : [""]
     nodepositions = @lift $(WF[:nodelabels]) ? coords : (dim < 3 ? Point2f[Point2f((0,0))] : Point3f[Point3f((0,0,0))])
-    celllabels = @lift $(WF[:celllabels]) ? ["$i" for i in 1:Ferrite.getncells(grid)] : [""]
-    cellpositions = @lift $(WF[:celllabels]) ? [midpoint(cell,coords) for cell in Ferrite.getcells(grid)] : (dim < 3 ? [Point2f((0,0))] : [Point3f((0,0,0))])
+    celllabels = @lift $(WF[:celllabels]) ? ["$i" for i in 1:ncells] : [""]
+    cellpositions = @lift $(WF[:celllabels]) ? [midpoint(cell,coords) for cell in allcells] : (dim < 3 ? [Point2f((0,0))] : [Point3f((0,0,0))])
     #cellsetsplot
-    if isconcretetype(grid.cells)
-        dh = Ferrite.DofHandler(grid)
-    else
-        dh = Ferrite.MixedDofHandler(grid)
-    end
+
+    dh = Ferrite.DofHandler(grid)
     cellsets = grid.cellsets
     cellset_to_value = Dict{String,Int}()
     for (cellsetidx,(cellsetname,cellset)) in enumerate(cellsets)
         cellset_to_value[cellsetname] = cellsetidx
     end
-    cellset_u = zeros(Ferrite.getncells(grid))
-    for (cellidx,cell) in enumerate(Ferrite.getcells(grid))
+    cellset_u = zeros(ncells)
+    for (cellidx,cell) in enumerate(allcells)
         for (cellsetname,cellsetvalue) in cellset_to_value
             if cellidx in cellsets[cellsetname]
                 cellset_u[cellidx] = cellsetvalue
@@ -283,10 +287,10 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:Ferrite.AbstractGrid{dim}}}) where
     plotter = MakiePlotter(dh,cellset_u)
     cellset_u =  reshape(transfer_scalar_celldata(plotter, cellset_u; process=identity), num_vertices(plotter))
     colorrange = isempty(cellset_to_value) ? (0,1) : (0,maximum(values(cellset_to_value)))
-    Makie.mesh!(WF, plotter.mesh, color=cellset_u, shading=false, scale_plot=false, colormap=:darktest, visible=WF[:cellsets])
-    Makie.text!(WF,nodelabels, position=nodepositions, textsize=WF[:textsize], offset=WF[:offset],color=WF[:nodelabelcolor])
-    Makie.text!(WF,celllabels, position=cellpositions, textsize=WF[:textsize], color=WF[:celllabelcolor], align=(:center,:center))
-    Makie.linesegments!(WF,lines,color=WF[:color], strokewidth=WF[:strokewidth], visible=WF[:visible])
+    Makie.mesh!(WF, plotter.mesh, color=cellset_u, shading=Makie.NoShading, colormap=:darktest, visible=WF[:cellsets])
+    Makie.text!(WF,nodelabels, position=nodepositions, fontsize=WF[:fontsize], offset=WF[:offset],color=WF[:nodelabelcolor])
+    Makie.text!(WF,celllabels, position=cellpositions, fontsize=WF[:fontsize], color=WF[:celllabelcolor], align=(:center,:center))
+    Makie.linesegments!(WF,lines,color=WF[:color], linewidth=WF[:linewidth], visible=WF[:visible])
 end
 
 """
@@ -299,8 +303,7 @@ values are transformed to a scalar based on `process` which defaults to the magn
 
 - `field = :default`
 - `process = postprocess`
-- `scale_plot = false`
-- `shading = false`
+- `shading = Makie.NoShading`
 - `colormap = :cividis`
 - `colorrange=Makie.automatic`
 - `nan_color::Union{Symbol, <:Colorant}=:red`
@@ -309,8 +312,7 @@ values are transformed to a scalar based on `process` which defaults to the magn
     Attributes(
     field = :default,
     process = postprocess,
-    scale_plot = false,
-    shading = false,
+    shading = Makie.NoShading,
     colormap = :cividis,
     colorrange=Makie.automatic,
     nan_color=:red,
@@ -330,7 +332,7 @@ function Makie.plot!(SF::Surface{<:Tuple{<:MakiePlotter{2}}})
     coords = @lift begin
         Point3f[Point3f(coord[1], coord[2], $(solution)[idx]) for (idx, coord) in enumerate(plotter.physical_coords)]
     end
-    return Makie.mesh!(SF, coords, plotter.vis_triangles, color=solution, scale_plot=SF[:scale_plot], shading=SF[:shading], colormap=SF[:colormap], colorrange=SF[:colorrange], nan_color=SF[:nan_color])
+    return Makie.mesh!(SF, coords, plotter.vis_triangles, color=solution, shading=SF[:shading], colormap=SF[:colormap], colorrange=SF[:colorrange], nan_color=SF[:nan_color])
 end
 
 """
@@ -361,26 +363,29 @@ the arrows are unicolored. Otherwise the color corresponds to the magnitude, or
     )
 end
 
-function Makie.plot!(AR::Arrows{<:Tuple{<:MakiePlotter{dim}}}) where dim
+function Makie.plot!(AR::Arrows{<:Tuple{<:MakiePlotter{sdim}}}) where sdim
     plotter = AR[1][]
     solution = @lift begin
         if $(AR[:field]) === :default
             field_name = Ferrite.getfieldnames(plotter.dh)[1]
-            @assert Ferrite.getfielddim(plotter.dh,field_name) > 1
+            field_dim = Ferrite.n_components(plotter.dh, field_name)
+            @assert field_dim == sdim "Dimension of field $field_name is $field_dim does not match spatial dimension $sdim"
             transfer_solution(plotter,$(plotter.u); field_name=field_name, process=identity)
         else
-            @assert Ferrite.getfielddim(plotter.dh,$(AR[:field])) > 1
+            field_name = $(AR[:field])
+            field_dim = Ferrite.n_components(plotter.dh, field_name)
+            @assert field_dim == sdim "Dimension of field $field_name is $field_dim does not match spatial dimension $sdim"
             transfer_solution(plotter,$(plotter.u); field_name=$(AR[:field]), process=identity)
         end
     end
-    if dim  == 2
+    if sdim  == 2
         ns = @lift([Vec2f(i) for i in eachrow($(solution))])
         lengths = @lift($(AR[:color])===:default ? $(AR[:process]).($(ns)) : ones(length($(ns)))*$(AR[:color]))
-    elseif dim  == 3
+    elseif sdim  == 3
         ns = @lift([Vec3f(i) for i in eachrow($(solution))])
         lengths = @lift($(AR[:color])===:default ? $(AR[:process]).($(ns)) : ones(length($(ns)))*$(AR[:color]))
     else
-        error("Arrows plots are only available in dim ≥ 2")
+        error("Arrows plots are only available in dimension ≥ 2")
     end
     Makie.arrows!(AR, plotter.physical_coords, ns, arrowsize=AR[:arrowsize], colormap=AR[:colormap], color=lengths, lengthscale=AR[:lengthscale])
 end
@@ -392,10 +397,10 @@ end
     elementinfo(cell::Type{AbstractCell}; kwargs...)
 
 - `plotnodes=true` controls if nodes of element are plotted
-- `strokewidth=2` strokwidth of faces/edges
+- `linewidth=2` strokwidth of faces/edges
 - `color=theme(scene, :linecolor)`
 - `markersize=30` size of the nodes
-- `textsize=60` textsize of node-, edges- and facelabels
+- `fontsize=60` fontsize of node-, edges- and facelabels
 - `nodelabels=true` switch that controls plotting of nodelabels
 - `nodelabelcolor=:darkred`
 - `nodelabeloffset=(0.0,0.0)` offset of the nodelabel text relative to its associated node
@@ -410,13 +415,16 @@ end
 @recipe(Elementinfo) do scene
     Attributes(
     plotnodes=true,
-    strokewidth=theme(scene, :linewidth),
+    linewidth=theme(scene, :linewidth),
     color=theme(scene, :linecolor),
     markersize=theme(scene, :markersize),
-    textsize=60,
+    fontsize=60,
+    vertexlabels=true,
+    vertexlabelcolor=:darkred,
+    vertexlabeloffset=(0.0,0.0),
     nodelabels=true,
     nodelabelcolor=:darkred,
-    nodelabeloffset=(0.0,0.0),
+    nodelabeloffset=(0.0,20.0),
     facelabels=true,
     facelabelcolor=:darkgreen,
     facelabeloffset=(-40,0),
@@ -427,61 +435,136 @@ end
     )
 end
 
-function Makie.plot!(Ele::Elementinfo{<:Tuple{<:Ferrite.Interpolation{dim,refshape}}}) where {dim,refshape}
-    ip = Ele[1][]
-    elenodes = Ferrite.reference_coordinates(ip) |> x->reshape(reinterpret(Float64,x),(dim,length(x)))'
+function Makie.plot!(Ele::Elementinfo{<:Tuple{<:Ferrite.AbstractCell{refshape}}}) where {refshape}
+    cell = Ele[1][]
+    dim = Ferrite.getrefdim(cell)
+
+    # Draw element outline
+    gip = Ferrite.geometric_interpolation(cell)
+    elenodes = Ferrite.reference_coordinates(gip) |> x->reshape(reinterpret(Float64,x),(dim,length(x)))'
     dim > 2 ? (lines = Point3f[]) : (lines = Point2f[])
-    facenodes = Ferrite.faces(ip)
-    if dim == 2
-        append!(lines, [elenodes[e,:] for boundary in facenodes for e in boundary[1:2]]) # 1:2 because higher order node in the middle
-    else
-        edgenodes = Ferrite.edges(ip)
-        order = Ferrite.getorder(ip)
-        #TODO remove the index monstrosity below after edges are defined consistently see https://github.com/Ferrite-FEM/Ferrite.jl/issues/520
-        append!(lines, [elenodes[e,:] for boundary in edgenodes for e in boundary[1:((refshape == Ferrite.RefCube) ? 1 : (order > 1 ? 2 : 1)):((refshape == Ferrite.RefCube) ? 2 : end)]]) # 1:2 because higher order node in the middle
-    end
-    boundaryentities = dim == 2 ? facenodes : edgenodes
-    #plot element boundary
-    Makie.linesegments!(Ele,lines,color=Ele[:color], linewidth=Ele[:strokewidth])
-    for (id,face) in enumerate(facenodes)
-        idx = 0
-        if refshape == Ferrite.RefCube && dim == 3
-            idx = 4
-        elseif refshape == Ferrite.RefTetrahedron && dim == 3
-            idx = 3
-        else
-            idx = 2
+    for edgenodes in Ferrite.edgedof_indices(gip)
+        append!(lines, [elenodes[edgenodes[1],:], elenodes[edgenodes[2],:]]) # The convention in Ferrite is that the first two edge nodes are always associated to the vertices
+    end
+    Makie.linesegments!(Ele,lines,color=Ele[:color], linewidth=Ele[:linewidth])
+
+    # Draw its nodes
+    Makie.scatter!(Ele,elenodes,markersize=Ele[:markersize], color=Ele[:color], visible=Ele[:plotnodes])
+    nodelabels = @lift $(Ele[:nodelabels]) ? ["N$i" for i in 1:size(elenodes,1)] : [""]
+    nodepositions = @lift $(Ele[:nodelabels]) ? [dim < 3 ? Point2f(row) : Point3f(row) for row in eachrow(elenodes)] : (dim < 3 ? [Point2f((0,0))] : [Point3f((0,0,0))])
+    Makie.text!(Ele,nodelabels, position=nodepositions, fontsize=Ele[:fontsize], offset=Ele[:nodelabeloffset],color=Ele[:nodelabelcolor],font=Ele[:font])
+
+    # Annotate element vertices
+    if dim ≥ 1 && Ele[:vertexlabels][]
+        for (id,vertexnodes) in enumerate(Ferrite.vertexdof_indices(gip))
+            position = if dim == 3
+                Point3f(elenodes[vertexnodes[1],:])
+            elseif dim == 2
+                Point2f(elenodes[vertexnodes[1],:])
+            end
+            Makie.text!(Ele,"V$id", position=position, fontsize=Ele[:fontsize], offset=Ele[:edgelabeloffset],color=Ele[:edgelabelcolor],visible=Ele[:edgelabels],font=Ele[:font])
         end
-        position = zeros(dim)
-        for i in 1:idx
-            position += elenodes[face[i],:]
+    end
+
+    # Annotate element edges
+    if dim ≥ 2 && Ele[:edgelabels][]
+        for (id,edgenodes) in enumerate(Ferrite.edgedof_indices(gip))
+            position = if dim == 3
+                Point3f((elenodes[edgenodes[1],:] + elenodes[edgenodes[2],:])*0.5)
+            elseif dim == 2
+                Point2f((elenodes[edgenodes[1],:] + elenodes[edgenodes[2],:])*0.5)
+            end
+            Makie.text!(Ele,"E$id", position=position, fontsize=Ele[:fontsize], offset=Ele[:edgelabeloffset],color=Ele[:edgelabelcolor],visible=Ele[:edgelabels],font=Ele[:font])
         end
-        position ./= idx
-        position = dim == 2 ? Point2f(position) : Point3f(position)
-        Makie.text!(Ele,"$id", position=position, textsize=Ele[:textsize], offset=Ele[:facelabeloffset],color=Ele[:facelabelcolor],visible=Ele[:facelabels],font=Ele[:font])
-    end
-    if dim == 3
-        for (id,edge) in enumerate(edgenodes)
-            position = Point3f((elenodes[edge[1],:] + elenodes[refshape==Ferrite.RefCube ? edge[2] : edge[end],:])*0.5)
-            t = Makie.text!(Ele,"$id", position=position, textsize=Ele[:textsize], offset=Ele[:edgelabeloffset],color=Ele[:edgelabelcolor],visible=Ele[:edgelabels],align=(:center,:center),font=Ele[:font])
-            # Boundingbox can't switch currently from pixelspace to "coordinate" space in recipes
-            #bb = Makie.boundingbox(t)
-            #Makie.wireframe!(Ele,bb,space=:pixel)
+    end
+
+    # Annotate element faces
+    if dim ≥ 3 && Ele[:facelabels][]
+        for face_index ∈ 1:Ferrite.nfaces(cell)
+            linear_face = linear_face_cell(cell, face_index)
+            gip_face = Ferrite.geometric_interpolation(linear_face)
+            vertexdof_list = Ferrite.vertexdof_indices(gip_face)
+            refpositions   = Ferrite.reference_coordinates(gip_face)
+            facenodes      = [Ferrite.facet_to_element_transformation(refposition, refshape, face_index) for refposition in refpositions] |> x->reshape(reinterpret(Float64,x),(dim,length(x)))'        
+            position = zeros(dim)
+            for i in 1:length(vertexdof_list)
+                position += facenodes[vertexdof_list[i][1],:]
+            end
+            position ./= length(vertexdof_list)
+            position = dim == 2 ? Point2f(position) : Point3f(position)
+            Makie.text!(Ele,"F$face_index", position=position, fontsize=Ele[:fontsize], offset=Ele[:facelabeloffset],color=Ele[:facelabelcolor],visible=Ele[:facelabels],font=Ele[:font])
         end
     end
-    #plot the nodes
+end
+
+function Makie.plot!(Ele::Elementinfo{<:Tuple{<:Ferrite.Interpolation{refshape}}}) where {refshape}
+    ip = Ele[1][]
+    dim = Ferrite.getrefdim(ip)
+
+    # Draw element boundary
+    gip = Ferrite.default_geometric_interpolation(ip)
+    geonodes = Ferrite.reference_coordinates(gip) |> x->reshape(reinterpret(Float64,x),(dim,length(x)))'
+    dim > 2 ? (lines = Point3f[]) : (lines = Point2f[])
+    for edgenodes in Ferrite.edgedof_indices(gip.ip)
+        append!(lines, [geonodes[edgenodes[1],:], geonodes[edgenodes[2],:]]) 
+    end
+    Makie.linesegments!(Ele,lines,color=Ele[:color], linewidth=Ele[:linewidth])
+
+    # Draw the interpolation nodes and its indices
+    elenodes = Ferrite.reference_coordinates(ip) |> x->reshape(reinterpret(Float64,x),(dim,length(x)))'
     Makie.scatter!(Ele,elenodes,markersize=Ele[:markersize], color=Ele[:color], visible=Ele[:plotnodes])
-    #set up nodelabels
-    nodelabels = @lift $(Ele[:nodelabels]) ? ["$i" for i in 1:size(elenodes,1)] : [""]
+    nodelabels = @lift $(Ele[:nodelabels]) ? ["D$i" for i in 1:size(elenodes,1)] : [""]
     nodepositions = @lift $(Ele[:nodelabels]) ? [dim < 3 ? Point2f(row) : Point3f(row) for row in eachrow(elenodes)] : (dim < 3 ? [Point2f((0,0))] : [Point3f((0,0,0))])
-    #set up celllabels
-    Makie.text!(Ele,nodelabels, position=nodepositions, textsize=Ele[:textsize], offset=Ele[:nodelabeloffset],color=Ele[:nodelabelcolor],font=Ele[:font])
-    #plot edges (3D) /faces (2D) of the mesh
-    Makie.linesegments!(Ele,lines,color=Ele[:color], linewidth=Ele[:strokewidth])
+    Makie.text!(Ele,nodelabels, position=nodepositions, fontsize=Ele[:fontsize], offset=Ele[:nodelabeloffset],color=Ele[:nodelabelcolor],font=Ele[:font])
+
+    # Annotate element vertices
+    if dim ≥ 1 && Ele[:vertexlabels][]
+        for (id,vertexnodes) in enumerate(Ferrite.vertexdof_indices(gip.ip))
+            position = if dim == 3
+                Point3f(elenodes[vertexnodes[1],:])
+            elseif dim == 2
+                Point2f(elenodes[vertexnodes[1],:])
+            end
+            Makie.text!(Ele,"V$id", position=position, fontsize=Ele[:fontsize], offset=Ele[:edgelabeloffset],color=Ele[:edgelabelcolor],visible=Ele[:edgelabels],font=Ele[:font])
+        end
+    end
+    # Annotate element edges
+    if dim ≥ 2 && Ele[:edgelabels][]
+        for (id,edgenodes) in enumerate(Ferrite.edgedof_indices(gip.ip))
+            position = if dim == 3
+                Point3f((elenodes[edgenodes[1],:] + elenodes[edgenodes[2],:])*0.5)
+            elseif dim == 2
+                Point2f((elenodes[edgenodes[1],:] + elenodes[edgenodes[2],:])*0.5)
+            end
+            Makie.text!(Ele,"E$id", position=position, fontsize=Ele[:fontsize], offset=Ele[:edgelabeloffset],color=Ele[:edgelabelcolor],visible=Ele[:edgelabels],font=Ele[:font])
+        end
+    end
+
+    # Annotate element faces
+    if dim ≥ 3 && Ele[:facelabels][]
+        for face_index ∈ 1:Ferrite.nfaces(ip)
+            linear_face = linear_face_cell(refshape, face_index)
+            gip_face = Lagrange{linear_face,1}()
+            vertexdof_list = Ferrite.vertexdof_indices(gip_face)
+            refpositions   = Ferrite.reference_coordinates(gip_face)
+            facenodes      = [Ferrite.facet_to_element_transformation(refposition, refshape, face_index) for refposition in refpositions] |> x->reshape(reinterpret(Float64,x),(dim,length(x)))'        
+            position = zeros(dim)
+            for i in 1:length(vertexdof_list)
+                position += facenodes[vertexdof_list[i][1],:]
+            end
+            position ./= length(vertexdof_list)
+            position = dim == 2 ? Point2f(position) : Point3f(position)
+            Makie.text!(Ele,"F$face_index", position=position, fontsize=Ele[:fontsize], offset=Ele[:facelabeloffset],color=Ele[:facelabelcolor],visible=Ele[:facelabels],font=Ele[:font])
+        end
+    end
 end
 
-Makie.convert_arguments(P::Type{<:Elementinfo}, cell::C) where C<:Ferrite.AbstractCell = (Ferrite.default_interpolation(typeof(cell)),)
-Makie.convert_arguments(P::Type{<:Elementinfo}, celltype::Type{C}) where C<:Ferrite.AbstractCell = (Ferrite.default_interpolation(celltype),)
+function Makie.convert_arguments(P::Type{<:Elementinfo}, celltype::Type{C}) where C<:Ferrite.AbstractCell 
+    gip = geometric_interpolation(C)
+    nnodes = getnbasefunctions(gip)
+    nodes  = ntuple(x->1,nnodes)
+    return (celltype(nodes),)
+end
 Makie.convert_arguments(P::Type{<:Elementinfo}, iptype::Type{IP}) where IP<:Ferrite.Interpolation = (iptype(),)
 
 """
@@ -504,12 +587,12 @@ function ferriteviewer(plotter::MakiePlotter{dim}) where dim
 
     #setting up various sliders
     markerslider = Slider(fig, range = 0:1:100, startvalue=5)
-    strokewidthslider = Slider(fig, range = 0:1:10, startvalue=1)
+    linewidthslider = Slider(fig, range = 0:1:10, startvalue=1)
     markersize = lift(x->x,markerslider.value)
-    strokewidth = lift(x->x,strokewidthslider.value)
+    linewidth = lift(x->x,linewidthslider.value)
 
     #plot the fe-mesh
-    wireframep = wireframe!(plotter,markersize=markersize,strokewidth=strokewidth,deformation_field= @lift $(toggles[2].active) ? $(deformation_field) : :default)
+    wireframep = wireframe!(plotter,markersize=markersize,linewidth=linewidth,deformation_field= @lift $(toggles[2].active) ? $(deformation_field) : :default)
     #connect fe-mesh plot to the toggle
     connect!(wireframep.visible,toggles[1].active)
     connect!(wireframep.nodelabels,toggles[3].active)
@@ -523,7 +606,7 @@ function ferriteviewer(plotter::MakiePlotter{dim}) where dim
     #align all menus as a vgrid under each other
     fig[1,3] = vgrid!(grid!(hcat(toggles,labels), tellheight=false),
                       Label(fig,"nodesize",width=nothing), markerslider,
-                      Label(fig,"strokewidth",width=nothing), strokewidthslider,
+                      Label(fig,"linewidth",width=nothing), linewidthslider,
                       Label(fig,"processing function",width=nothing), menu_process,
                       Label(fig,"field",width=nothing), menu_field,
                       Label(fig, "deformation field",width=nothing),menu_deformation_field,
@@ -562,7 +645,25 @@ end
 
 ####### One Shot Methods #######
 const FerriteVizPlots = Union{Type{<:Wireframe},Type{<:SolutionPlot},Type{<:Arrows},Type{<:Surface}}
+# We default with our axis choice to the spatial dimension of the problem
+function Makie.args_preferred_axis(a, b::Union{MakiePlotter{sdim},Grid{sdim}}, args...) where {sdim}
+    if sdim ≤ 2
+        return Makie.Axis
+    else
+        return Makie.LScene
+    end
+end
+function Makie.args_preferred_axis(a::Type{<:Elementinfo}, ip_or_cell)
+    dim = Ferrite.getrefdim(ip_or_cell)
+    if dim ≤ 2
+        return Makie.Axis
+    else
+        return Makie.LScene
+    end
+end
+# Surface plots are special, as they are 2D problems which are deformed into the third dimension
+Makie.args_preferred_axis(a::Type{<:Surface}, b::Union{MakiePlotter{sdim},Grid{sdim}}, args...) where {sdim} = Makie.LScene
 
-function Makie.convert_arguments(P::FerriteVizPlots, dh::Ferrite.AbstractDofHandler, u::Vector)
+function Makie.convert_arguments(P::FerriteVizPlots, dh::Ferrite.AbstractDofHandler, u::AbstractVector)
     return (MakiePlotter(dh,u),)
 end
diff --git a/src/utils.jl b/src/utils.jl
index ad99d91..b342174 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,21 +1,12 @@
-# Helper... Refactoring needed.
-function getfieldinterpolation(dh::Ferrite.DofHandler, field_name::Symbol)
-    field_idx = indexin(dh.field_names, [:b])
-    field_idx == nothing && error("did not find field $field_name")
-    dh.field_interpolations[field_idx]
-end
-
-function field_offset(dh::Ferrite.DofHandler, field_idx::Int)
-    offset = 0
-    for i in 1:(field_idx-1)
-        offset += Ferrite.getnbasefunctions(dh.field_interpolations[i])::Int * dh.field_dims[i]
-    end
-    return offset
-end
+# Helper for shape of geometry
+const TriangularCell = Ferrite.AbstractCell{<:Ferrite.RefTriangle}
+const QuadrilateralCell = Ferrite.AbstractCell{<:Ferrite.RefQuadrilateral}
+const HexahedralCell = Ferrite.AbstractCell{<:Ferrite.RefHexahedron}
+const TetrahedralCell = Ferrite.AbstractCell{<:Ferrite.RefTetrahedron}
 
-Ferrite.vertices(cell::Ferrite.Cell{3,3,1}) = cell.nodes
-
-Ferrite.default_interpolation(::Type{Ferrite.Cell{3,3,1}}) = Ferrite.Lagrange{2,Ferrite.RefTetrahedron,1}()
+# FIXME this is not correct in general and may lead to loss of details
+get_gradient_interpolation(::Lagrange{shape,order}) where {sdim,shape<:Ferrite.AbstractRefShape{sdim},order} = VectorizedInterpolation{sdim}(DiscontinuousLagrange{shape,order-1}())
+get_gradient_interpolation_type(::Type{Lagrange{shape,order}}) where {sdim,shape<:Ferrite.AbstractRefShape{sdim},order} = VectorizedInterpolation{sdim,shape,order-1,DiscontinuousLagrange{shape,order-1}}
 
 """
     linear_face_cell(cell::Ferrite.Cell, local_face_idx::Int)
@@ -25,11 +16,10 @@ Get the geometrically linear face of a given cell.
 !!! warning
     This may extracts the face spanned by the vertices, not the actual face!
 """
-linear_face_cell(cell::Ferrite.Cell{3,N,4}, local_face_idx::Int) where N =  Ferrite.Cell{3,3,1}(Ferrite.faces(cell)[local_face_idx])
-linear_face_cell(cell::Ferrite.Cell{3,N,6}, local_face_idx::Int) where N = Ferrite.Quadrilateral3D(Ferrite.faces(cell)[local_face_idx])
-
-# Obtain the face interpolation on regular geometries.
-getfaceip(ip::Ferrite.Interpolation{dim, shape, order}, local_face_idx::Int) where {dim, shape <: Union{Ferrite.RefTetrahedron, Ferrite.RefCube}, order} = Ferrite.getlowerdim(ip)
+linear_face_cell(cell::TetrahedralCell, local_face_idx::Int) = Triangle(Ferrite.faces(cell)[local_face_idx])
+linear_face_cell(cell::HexahedralCell, local_face_idx::Int)  = Quadrilateral(Ferrite.faces(cell)[local_face_idx])
+linear_face_cell(cell::Type{RefHexahedron}, local_face_idx::Int)  = RefQuadrilateral
+linear_face_cell(cell::Type{RefTetrahedron}, local_face_idx::Int) = RefTriangle
 
 struct MakiePlotter{dim,DH<:Ferrite.AbstractDofHandler,T1,TOP<:Union{Nothing,Ferrite.AbstractTopology},T2,M,TRI} <: AbstractPlotter
     dh::DH
@@ -62,10 +52,10 @@ The triangulation acts as a "L2" triangulation, i.e. the nodes which are shared
 """
 function MakiePlotter(dh::Ferrite.AbstractDofHandler, u::Vector, topology::TOP) where {TOP<:Union{Nothing,Ferrite.AbstractTopology}}
     cells = Ferrite.getcells(dh.grid)
-    dim = Ferrite.getdim(dh.grid)
+    dim = Ferrite.getspatialdim(dh.grid)
     visible = zeros(Bool,length(cells))
     if dim > 2
-        boundaryfaces = findall(isempty,topology.face_neighbor)
+        boundaryfaces = findall(isempty,topology.face_face_neighbor)
         boundaryelements = Ferrite.getindex.(boundaryfaces,1)
     else
         boundaryelements = collect(1:Ferrite.getncells(dh.grid))
@@ -85,7 +75,7 @@ function MakiePlotter(dh::Ferrite.AbstractDofHandler, u::Vector, topology::TOP)
     triangles = Matrix{Int}(undef, num_triangles, 3)
     triangle_cell_map = Vector{Int}(undef, num_triangles)
     physical_coords = Vector{GeometryBasics.Point{dim,Float32}}(undef, num_triangles*3)
-    gridnodes = [GeometryBasics.Point{dim,Float32}(node.data) for node in Ferrite.getcoordinates.(Ferrite.getnodes(dh.grid))]
+    gridnodes = [GeometryBasics.Point{dim,Float32}(node.data) for node in Ferrite.get_node_coordinate.(Ferrite.getnodes(dh.grid))]
     reference_coords = Matrix{Float64}(undef, num_triangles*3,dim)
 
     # Decompose does the heavy lifting for us
@@ -106,7 +96,7 @@ function MakiePlotter(dh::Ferrite.AbstractDofHandler, u::Vector, topology::TOP)
     mesh = GeometryBasics.Mesh(physical_coords_m,vis_triangles)
     return MakiePlotter{dim,typeof(dh),eltype(u),typeof(topology),Float32,typeof(mesh),eltype(vis_triangles)}(dh,Observable(u),topology,visible,gridnodes,physical_coords,physical_coords_m,all_triangles,vis_triangles,triangle_cell_map,cell_triangle_offsets,reference_coords,mesh)
 end
-MakiePlotter(dh,u) = MakiePlotter(dh,u,Ferrite.getdim(dh.grid) > 2 ? Ferrite.ExclusiveTopology(dh.grid.cells) : nothing)
+MakiePlotter(dh,u) = MakiePlotter(dh,u,Ferrite.getspatialdim(dh.grid) > 2 ? Ferrite.ExclusiveTopology(dh.grid) : nothing)
 
 """
     ClipPlane{T}(normal, distance_to_origin)
@@ -123,9 +113,9 @@ struct ClipPlane{T}
 end
 
 # Binary decision function to clip a cell with a plane for the crinkle clip.
-function (plane::ClipPlane)(grid, cellid)
-    cell = grid.cells[cellid]
-    coords = Ferrite.getcoordinates.(Ferrite.getnodes(grid)[[cell.nodes...]])
+function (plane::ClipPlane)(grid::Ferrite.AbstractGrid, cellid::Int)
+    cell = getcells(grid, cellid)
+    coords = Ferrite.get_node_coordinate.(Ferrite.getnodes(grid)[[cell.nodes...]])
     for coord ∈ coords
         if coord ⋅ plane.normal > plane.distance
             return false
@@ -197,13 +187,13 @@ Total number of vertices
 """
 num_vertices(p::MakiePlotter) = size(p.physical_coords,1)
 
-# TODO this looks faulty...think harder.
-# Helper to count triangles e.g. for preallocations.
-ntriangles(cell::Ferrite.AbstractCell{2,N,3}) where {N} = 1 # Tris in 2D
-ntriangles(cell::Ferrite.AbstractCell{3,3,1}) = 1 # Tris in 3D
-ntriangles(cell::Ferrite.AbstractCell{dim,N,4}) where {dim,N} = 4 # Quads in 2D and 3D
-ntriangles(cell::Ferrite.AbstractCell{3,N,1}) where N = 4 # Tets as a special case of a Quad, obviously :)
-ntriangles(cell::Ferrite.AbstractCell{3,N,6}) where N = 6*4 # Hex
+"""
+Helper to count triangles in static triangulations, e.g. for preallocations.
+"""
+ntriangles(cell::TriangularCell) = 1
+ntriangles(cell::QuadrilateralCell) = 4
+ntriangles(cell::TetrahedralCell)  = 4
+ntriangles(cell::HexahedralCell) = 6*4
 
 """
 Get the vertices represented as a list of coordinates of a cell.
@@ -211,7 +201,7 @@ Get the vertices represented as a list of coordinates of a cell.
 !!! details
     **TODO** refactor into Ferrite core.
 """
-function vertices(grid::Ferrite.AbstractGrid, cell::Ferrite.AbstractCell{dim,N,M}) where {dim,N,M}
+function vertices(grid::Ferrite.AbstractGrid, cell::Ferrite.AbstractCell)
     Ferrite.getnodes(grid)[[Ferrite.vertices(cell)...]]
 end
 
@@ -220,10 +210,12 @@ end
 
 Decompose a triangle into a coordinates and a triangle index list to disconnect it properly. Guarantees to preserve orderings and orientations.
 """
-function decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Union{Ferrite.AbstractCell{2,N,3}, Ferrite.AbstractCell{3,3,1}}) where {N}
+function decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::TriangularCell)
+    ip_geo = Ferrite.geometric_interpolation(typeof(cell))
+    refcoords = Ferrite.reference_coordinates(ip_geo)
     for (i,v) in enumerate(vertices(grid, cell))
-        coord_matrix[coord_offset] = GeometryBasics.Point(Ferrite.getcoordinates(v)...)
-        ref_coord_matrix[coord_offset,1:2] = Ferrite.reference_coordinates(Ferrite.default_interpolation(typeof(cell)))[i]
+        coord_matrix[coord_offset] = GeometryBasics.Point(Ferrite.get_node_coordinate(v)...)
+        ref_coord_matrix[coord_offset,1:2] = refcoords[i]
         triangle_matrix[triangle_offset,i] = coord_offset
         coord_offset+=1
     end
@@ -259,19 +251,22 @@ Decompose a quadrilateral into a coordinates and a triangle index list to discon
     ```
     where A=(1,2,5),B=(2,3,5),C=(3,4,5),D=(4,1,5) are the generated triangles in this order.
 """
-function decompose!(coord_offset, coord_matrix::Vector{Point{space_dim,T}}, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Union{Ferrite.AbstractCell{2,N,4}, Ferrite.AbstractCell{3,4,1}}) where {N,space_dim,T}
+function decompose!(coord_offset, coord_matrix::Vector{Point{sdim,T}}, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::QuadrilateralCell) where {sdim, T}
     # A bit more complicated. The default diagonal decomposition into 2 triangles is missing a solution mode.
     # To resolve this we make a more expensive decomposition in 4 triangles which correctly displays the solution mode (in linear problems).
     coord_offset_initial = coord_offset
     vts = vertices(grid, cell)
 
     # Compute coordinate of vertex 5
-    center = zeros(space_dim)
+    center = zeros(sdim)
     for v in vts
-        center += Ferrite.getcoordinates(v)
+        center += Ferrite.get_node_coordinate(v)
     end
     center /= 4.0
 
+    ip_geo = Ferrite.geometric_interpolation(typeof(cell))
+    refcoords = Ferrite.reference_coordinates(ip_geo)
+
     # Generate triangles in order
     for i = 1:length(vts)
         v1 = vts[i]
@@ -283,13 +278,13 @@ function decompose!(coord_offset, coord_matrix::Vector{Point{space_dim,T}}, ref_
         v2 = vts[i2]
 
         # current vertex
-        coord_matrix[coord_offset] = GeometryBasics.Point(Ferrite.getcoordinates(v1)...)
-        ref_coord_matrix[coord_offset, 1:2] = Ferrite.reference_coordinates(Ferrite.default_interpolation(typeof(cell)))[i]
+        coord_matrix[coord_offset] = GeometryBasics.Point(Ferrite.get_node_coordinate(v1)...)
+        ref_coord_matrix[coord_offset, 1:2] = refcoords[i]
         coord_offset+=1
 
         # next vertex in chain
-        coord_matrix[coord_offset] = GeometryBasics.Point(Ferrite.getcoordinates(v2)...)
-        ref_coord_matrix[coord_offset, 1:2] = Ferrite.reference_coordinates(Ferrite.default_interpolation(typeof(cell)))[i2]
+        coord_matrix[coord_offset] = GeometryBasics.Point(Ferrite.get_node_coordinate(v2)...)
+        ref_coord_matrix[coord_offset, 1:2] = refcoords[i2]
         coord_offset+=1
 
         # center vertex (5)
@@ -306,36 +301,33 @@ function decompose!(coord_offset, coord_matrix::Vector{Point{space_dim,T}}, ref_
 end
 
 """
-    decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Ferrite.AbstractCell{3,N,M})
+    decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Ferrite.AbstractCell{<:Ferrite.AbstractRefShape{3}})
 
 Decompose volumetric objects via their faces.
 """
-function decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Ferrite.AbstractCell{3,N,M}) where {N,M}
-    # Just 6 quadrilaterals :)
-    for face_index ∈ 1:M
+function decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Ferrite.AbstractCell{CellType}) where CellType<:Ferrite.AbstractRefShape{3}
+    for face_index ∈ 1:Ferrite.nfaces(cell)
         face_coord_offset = coord_offset
         (coord_offset, triangle_offset) = decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, linear_face_cell(cell, face_index))
         for ci ∈ face_coord_offset:(coord_offset-1)
-            new_coord = transfer_quadrature_face_to_cell(ref_coord_matrix[ci, 1:2], cell, face_index)
+            new_coord = Ferrite.facet_to_element_transformation(Tensors.Vec{2}(ref_coord_matrix[ci, 1:2]), CellType, face_index)
             ref_coord_matrix[ci, :] = new_coord
         end
     end
     (coord_offset, triangle_offset)
 end
 
-
-refshape(cell::Ferrite.AbstractCell) = typeof(Ferrite.default_interpolation(typeof(cell))).parameters[2]
-
 x₁(x) = x[1]
 x₂(x) = x[2]
 x₃(x) = x[3]
 l2(x) = LinearAlgebra.norm(x,2)
 l1(x) = LinearAlgebra.norm(x,1)
 
-midpoint(cell::Ferrite.AbstractCell{2,N,3}, points) where N = Point2f((1/3) * (points[cell.nodes[1]] + points[cell.nodes[2]] + points[cell.nodes[3]]))
-midpoint(cell::Ferrite.AbstractCell{2,N,4}, points) where N = Point2f(0.5 * (points[cell.nodes[1]] + points[cell.nodes[3]]))
-midpoint(cell::Ferrite.AbstractCell{3,N,4}, points) where N = Point3f((1/4) * (points[cell.nodes[1]] + points[cell.nodes[2]] + points[cell.nodes[3]] + points[cell.nodes[4]]))
-midpoint(cell::Ferrite.AbstractCell{3,N,6}, points) where N = Point3f(0.5 * (points[cell.nodes[1]] + points[cell.nodes[7]]))
+# TODO deal with nonlinear geometries
+midpoint(cell::TriangularCell, points) = Point2f((1/3) * (points[cell.nodes[1]] + points[cell.nodes[2]] + points[cell.nodes[3]]))
+midpoint(cell::QuadrilateralCell, points) = Point2f(0.5 * (points[cell.nodes[1]] + points[cell.nodes[3]]))
+midpoint(cell::TetrahedralCell, points) = Point3f((1/4) * (points[cell.nodes[1]] + points[cell.nodes[2]] + points[cell.nodes[3]] + points[cell.nodes[4]]))
+midpoint(cell::HexahedralCell, points) = Point3f(0.5 * (points[cell.nodes[1]] + points[cell.nodes[7]]))
 
 """
     postprocess(node_values::Vector{T}) -> T
@@ -351,25 +343,15 @@ function postprocess(node_values)
     end
 end
 
-function getfieldhandlers(dh::Ferrite.DofHandler,field_name)
-    names = Ferrite.getfieldnames(dh)
-    field_idx = Ferrite.find_field.((dh,),names)
-    ip_field = Ferrite.getfieldinterpolation.((dh,),field_idx)
-    field_dim_ = Ferrite.getfielddim.((dh,),field_idx)
-    return [Ferrite.FieldHandler([Ferrite.Field(fname,fip,fdim) for (fname,fip,fdim) in zip(names,ip_field,field_dim_)],Set(1:Ferrite.getncells(dh.grid)))]
-end
-
-function getfieldhandlers(dh::Ferrite.MixedDofHandler,field_name)
-    fhs = Ferrite.FieldHandler[]
-    for fh in dh.fieldhandlers
-        for field in fh.fields
-            if field.name == field_name
-                push!(fhs,fh)
-                break
-            end
+# TODO upstream to Ferrite.jl 
+function getsubdofhandlers(dh::Ferrite.DofHandler,field_name)
+    sdhs = SubDofHandler[]
+    for sdh in dh.subdofhandlers
+        if field_name ∈ sdh.field_names
+            push!(sdhs,sdh)
         end
     end
-    return fhs
+    return sdhs
 end
 
 """
@@ -384,43 +366,54 @@ function transfer_solution(plotter::MakiePlotter{dim,DH,T}, u::Vector; field_nam
     dh = plotter.dh
     grid = dh.grid
 
+    sdhs = getsubdofhandlers(dh,field_name)
+    sdh1 = first(sdhs)
+    ip_field = Ferrite.getfieldinterpolation(sdh1,field_name)
     # field related variables
-    field_dim = Ferrite.getfielddim(dh, field_name)
-    val_buffer = zeros(T,field_dim)
+    ref_shape = Ferrite.getrefshape(ip_field)
+    ref_dim = Ferrite.getrefdim(ref_shape)
+    # Create a coordinate inside the reference space to evaluate the shape function at
+    ξinside = Ferrite.Vec(ntuple(d->0.0,ref_dim))
+    val_buffer = Ferrite.reference_shape_value(ip_field, ξinside, 1)
+    # NOTE this does not work for fancy ansatz spaces where derivatives are mixed in (e.g. Hermite)
     val = process(val_buffer)
     _processreturn = length(process(val_buffer))
     data = fill(NaN, num_vertices(plotter),_processreturn)
-    for fh in getfieldhandlers(dh,field_name)
-        ip_field = Ferrite.getfieldinterpolation(fh,field_name)
-        cellset_ = collect(fh.cellset)
-        cell_geo_ref = Ferrite.getcells(grid, cellset_[1])
-        ntriangles(cell_geo_ref) == 0 && continue
-        ip_geo = Ferrite.default_interpolation(typeof(cell_geo_ref))
-        pv = Ferrite.PointScalarValues(ip_field, ip_geo)
-        _transfer_solution!(data,pv,fh,ip_geo,ip_field,cellset_,val_buffer,val,field_name,field_dim,plotter,u,process) #function barrier for ip_field and thus pointvalues
+
+    for sdh in getsubdofhandlers(dh,field_name)
+        cellset_ = collect(sdh.cellset)
+        ip_field = Ferrite.getfieldinterpolation(sdh,field_name)
+        ip_geo = Ferrite.geometric_interpolation(Ferrite.getcelltype(sdh))
+        pv = Ferrite.PointValues(ip_field, ip_geo; update_gradients=false)
+        field_dim = Ferrite.n_components(dh, field_name)
+        ref_shape = Ferrite.getrefshape(ip_field)
+        val_buffer = Ferrite.reference_shape_value(ip_field, ξinside, 1)
+        val = process(val_buffer)
+        _transfer_solution!(data,pv,sdh,ip_field,cellset_,val_buffer,val,field_name,field_dim,plotter,u,process) #function barrier for ip_field and thus pointvalues
     end
     return data
 end
 
-function _transfer_solution!(data,pv,fh,ip_geo,ip_field,cellset_,val_buffer,val,field_name,field_dim,plotter::MakiePlotter{dim,DH,T}, u::Vector, process::FUN) where {dim,DH<:Ferrite.AbstractDofHandler,T,FUN}
+function _transfer_solution!(data,pv,sdh,ip_field,cellset_,val_buffer,val,field_name,field_dim,plotter::MakiePlotter{dim,DH,T}, u::Vector, process::FUN) where {dim,DH<:Ferrite.AbstractDofHandler,T,FUN}
     n_vertices_per_tri = 3 # we have 3 vertices per triangle...
     dh = plotter.dh
     ref_coords = plotter.reference_coords
     grid = dh.grid
     # actual data
-    local_dof_range = Ferrite.dof_range(fh, field_name)
+    local_dof_range = Ferrite.dof_range(sdh, field_name)
     _processreturndim = length(process(val_buffer))
 
     cell_geo_ref = Ferrite.getcells(grid, cellset_[1])
 
-    Ferrite.reinit!(pv, Ferrite.getcoordinates(grid,cellset_[1]), Tensors.Vec{dim}(ref_coords[1,:]))
     n_basefuncs = Ferrite.getnbasefunctions(pv)
 
     _local_coords = Ferrite.getcoordinates(grid,cellset_[1])
     _local_celldofs = Ferrite.celldofs(dh,cellset_[1])
-    _celldofs_field = reshape(@view(_local_celldofs[local_dof_range]), (field_dim, n_basefuncs))
+    _celldofs_field = @view(_local_celldofs[local_dof_range])
     _local_ref_coords = Tensors.Vec{dim}(ref_coords[1,:])
 
+    Ferrite.reinit!(pv, _local_coords, Tensors.Vec{dim}(ref_coords[1,:]))
+
     # We just loop over all cells
     for (isvisible,(cell_idx,cell_geo)) in zip(plotter.visible,enumerate(Ferrite.getcells(dh.grid)))
         # Skip invisible cells and cells which are not in the current cellset
@@ -431,17 +424,14 @@ function _transfer_solution!(data,pv,fh,ip_geo,ip_field,cellset_,val_buffer,val,
         # Buffer cell data relevant for the current field to transfer
         Ferrite.getcoordinates!(_local_coords,grid,cell_idx)
         Ferrite.celldofs!(_local_celldofs,dh,cell_idx)
-        _celldofs_field = reshape(@view(_local_celldofs[local_dof_range]), (field_dim, n_basefuncs))
+        _celldofs_field = @view(_local_celldofs[local_dof_range])
 
         # Loop over the triangles of the cell and interpolate at the vertices
-        # TODO remove redundant function value calls
         for triangle_index in triangles_on_cell(plotter, cell_idx)
             for current_vertex_index in plotter.all_triangles[triangle_index]
                 _local_ref_coords = Tensors.Vec{dim}(@view(ref_coords[current_vertex_index,:]))
                 Ferrite.reinit!(pv, _local_coords, _local_ref_coords)
-                for d in 1:field_dim
-                    val_buffer[d] = Ferrite.function_value(pv, 1, @views(u[_celldofs_field[d,:]]))
-                end
+                val_buffer = Ferrite.function_value(pv, 1, @views(u[_celldofs_field]))
                 val = process(val_buffer)
                 for d in 1:_processreturndim
                     data[current_vertex_index, d] = val[d]
@@ -475,19 +465,6 @@ function transfer_scalar_celldata(plotter::MakiePlotter{dim,DH,T}, u::Vector; pr
     return data::Vector{T}
 end
 
-get_gradient_interpolation(::Ferrite.Lagrange{dim,shape,order}) where {dim,shape,order} = Ferrite.DiscontinuousLagrange{dim,shape,order-1}()
-get_gradient_interpolation_type(::Type{Ferrite.Lagrange{dim,shape,order}}) where {dim,shape,order} = Ferrite.DiscontinuousLagrange{dim,shape,order-1}
-# TODO remove if Knuth's PR on this gets merged (Ferrite PR 552)
-getgrid(dh::Ferrite.DofHandler) = dh.grid
-
-function ε(x::Vector{T}) where T
-    ngrad = length(x)
-    dim = isqrt(ngrad)
-    ∇u = Tensor{2,dim,T,ngrad}(x)
-    return symmetric(∇u)
-end
-
-
 """
     _tensorsjl_gradient_accessor(v::Tensors.Vec, field_dim_idx::Int, spatial_dim_idx::Int)
 
@@ -503,29 +480,30 @@ Compute the piecewise discontinuous gradient field for `field_name`. Returns the
 If the additional keyword argument `copy_fields` is provided with a non empty `Vector{Symbol}`, the corresponding fields of `dh` will be
 copied into the returned flux dof handler and flux dof value vector.
 """
-function interpolate_gradient_field(dh::Ferrite.DofHandler{spatial_dim}, u::AbstractVector, field_name::Symbol; copy_fields::Vector{Symbol}=Symbol[]) where {spatial_dim}
+function interpolate_gradient_field(dh::DofHandler{spatial_dim}, u::AbstractVector, field_name::Symbol; copy_fields::Vector{Symbol}=Symbol[]) where {spatial_dim}
+    @assert length(dh.subdofhandlers) == 1
+
     # Get some helpers
     field_idx = Ferrite.find_field(dh, field_name)
     ip = Ferrite.getfieldinterpolation(dh, field_idx)
+    field_dim = Ferrite.n_components(dh, field_idx)
 
     # Create dof handler for gradient field
-    dh_gradient = Ferrite.DofHandler(getgrid(dh))
+    dh_gradient = Ferrite.DofHandler(Ferrite.get_grid(dh))
     ip_gradient = get_gradient_interpolation(ip)
-    field_dim = Ferrite.getfielddim(dh,field_name)
-    Ferrite.add!(dh_gradient, :gradient, field_dim*spatial_dim, ip_gradient) # field dim × spatial dim components
+    add!(dh_gradient, :gradient, ip_gradient) # field dim × spatial dim components
     for fieldname in copy_fields
         _field_idx = Ferrite.find_field(dh, fieldname)
         _ip = Ferrite.getfieldinterpolation(dh, _field_idx)
-        _field_dim = Ferrite.getfielddim(dh,fieldname)
-        Ferrite.add!(dh_gradient, fieldname, _field_dim, _ip)
+        add!(dh_gradient, fieldname, _ip)
     end
     Ferrite.close!(dh_gradient)
 
     # FIXME this does not work for mixed grids
-    ip_geom = Ferrite.default_interpolation(typeof(Ferrite.getcells(getgrid(dh), 1)))
+    ip_geom = Ferrite.geometric_interpolation(typeof(Ferrite.getcells(Ferrite.get_grid(dh), 1)))
     ref_coords_gradient = Ferrite.reference_coordinates(ip_gradient)
-    qr_gradient = Ferrite.QuadratureRule{spatial_dim, refshape(Ferrite.getcells(getgrid(dh), 1)), Float64}(ones(length(ref_coords_gradient)), ref_coords_gradient)
-    cv = (field_dim == 1) ? Ferrite.CellScalarValues(qr_gradient, ip, ip_geom) : Ferrite.CellVectorValues(qr_gradient, ip, ip_geom)
+    qr_gradient = QuadratureRule{getrefshape(Ferrite.getcells(get_grid(dh), 1))}(ones(length(ref_coords_gradient)), ref_coords_gradient)
+    cv = CellValues(qr_gradient, ip, ip_geom)
 
     # Buffer for the dofs
     cell_dofs = zeros(Int, Ferrite.ndofs_per_cell(dh))
@@ -535,8 +513,9 @@ function interpolate_gradient_field(dh::Ferrite.DofHandler{spatial_dim}, u::Abst
     u_gradient = zeros(Ferrite.ndofs(dh_gradient))
     # In general uᵉ_gradient is an order 3 tensor [field_dim, spatial_dim, nqp]
     uᵉ_gradient = zeros(length(cell_dofs_gradient[Ferrite.dof_range(dh_gradient, :gradient)]))
-    uᵉshape = (spatial_dim, field_dim, Ferrite.getnquadpoints(cv))
+    uᵉshape = (spatial_dim, field_dim, getnquadpoints(qr_gradient))
     uᵉ_gradient_view = reshape(uᵉ_gradient, uᵉshape)
+    uᵉ = zeros(Ferrite.getnbasefunctions(ip))
 
     for (cell_num, cell) in enumerate(Ferrite.CellIterator(dh))
         # Get element dofs on parent field
@@ -548,7 +527,7 @@ function interpolate_gradient_field(dh::Ferrite.DofHandler{spatial_dim}, u::Abst
         Ferrite.reinit!(cv, cell)
 
         # Now we simply loop over all basis functions of the gradient field and evaluate the gradient
-        for i ∈ 1:Ferrite.getnquadpoints(cv)
+        for i ∈ 1:Ferrite.getnquadpoints(qr_gradient)
             uᵉgradi = Ferrite.function_gradient(cv, i, uᵉ)
             for ds in 1:spatial_dim
                 for df in 1:field_dim
@@ -571,9 +550,9 @@ end
 
 # maps the dof vector in nodal order, only needed for wireframe nodal deformation (since we display the original nodes)
 function dof_to_node(dh::Ferrite.AbstractDofHandler, u::Vector{T}; field_name=:u) where T
-    field_dim = Ferrite.getfielddim(dh, field_name)
+    field_dim = Ferrite.n_components(dh, field_name)
     data = fill(NaN, Ferrite.getnnodes(dh.grid), field_dim)
-    fhs = getfieldhandlers(dh,field_name)
+    fhs = getsubdofhandlers(dh,field_name)
 
     for fh in fhs
         dof_range_ = Ferrite.dof_range(fh, field_name)
@@ -590,34 +569,6 @@ function dof_to_node(dh::Ferrite.AbstractDofHandler, u::Vector{T}; field_name=:u
     return data
 end
 
-"""
-    transfer_quadrature_face_to_cell(point::AbstractVector, cell::Ferrite.AbstractCell{3,N,4}, face::Int)
-
-Mapping from 2D triangle to 3D face of a tetrahedon.
-"""
-function transfer_quadrature_face_to_cell(point::AbstractVector, cell::Ferrite.AbstractCell{3,N,4}, face::Int) where {N}
-    x,y = point
-    face == 1 && return [ 1-x-y,  y,  0]
-    face == 2 && return [ y,  0,  1-x-y]
-    face == 3 && return [ x,  y,  1-x-y]
-    face == 4 && return [ 0,  1-x-y,  y]
-end
-
-"""
-    transfer_quadrature_face_to_cell(point::AbstractVector, cell::Ferrite.AbstractCell{3,N,6}, face::Int) 
-
-Mapping from 2D quadrilateral to 3D face of a hexahedron.
-"""
-function transfer_quadrature_face_to_cell(point::AbstractVector, cell::Ferrite.AbstractCell{3,N,6}, face::Int) where {N}
-    x,y = point
-    face == 1 && return [ y,  x, -1]
-    face == 2 && return [ x, -1,  y]
-    face == 3 && return [ 1,  x,  y]
-    face == 4 && return [-x,  1,  y]
-    face == 5 && return [-1,  y,  x]
-    face == 6 && return [ x,  y,  1]
-end
-
 """
     uniform_refinement(plotter::MakiePlotter)
     uniform_refinement(plotter::MakiePlotter, num_refinements::Int)
@@ -725,3 +676,91 @@ function uniform_refinement(plotter::MakiePlotter{dim,DH,T1,TOP,T2,M,TRI}, num_r
     new_plotter = uniform_refinement(new_plotter, num_refinements-1)
     return new_plotter
 end
+
+# Foundation block for https://github.com/Ferrite-FEM/Ferrite.jl/issues/398 - Remove this below after the issue is resolved.
+##################################################
+# MatrixizedInterpolation{<:ScalarInterpolation} #
+##################################################
+abstract type MatrixInterpolation{vdim1, vdim2, refshape, order} <: Ferrite.Interpolation{refshape, order, Nothing} end
+
+struct MatrixizedInterpolation{vdim1, vdim2, refshape, order, SI <: ScalarInterpolation{refshape, order}} <: MatrixInterpolation{vdim1, vdim2, refshape,order}
+    ip::SI
+    function MatrixizedInterpolation{vdim1, vdim2}(ip::SI) where {vdim1, vdim2, refshape, order, SI <: Ferrite.ScalarInterpolation{refshape, order}}
+        return new{vdim1, vdim2, refshape, order, SI}(ip)
+    end
+end
+
+Ferrite.mapping_type(::MatrixizedInterpolation) = Ferrite.IdentityMapping()
+
+Ferrite.typeof_N(   ::Type{T}, ::MatrixizedInterpolation{dim,dim}, ::VectorizedInterpolation{dim, <: Ferrite.AbstractRefShape{dim}}) where {T, dim} = Tensor{2, dim, T}
+Ferrite.typeof_dNdx(::Type{T}, ::MatrixizedInterpolation{dim,dim}, ::VectorizedInterpolation{dim, <: Ferrite.AbstractRefShape{dim}}) where {T, dim} = Tensor{3, dim, T}
+Ferrite.typeof_dNdξ(::Type{T}, ::MatrixizedInterpolation{dim,dim}, ::VectorizedInterpolation{dim, <: Ferrite.AbstractRefShape{dim}}) where {T, dim} = Tensor{3, dim, T}
+
+Ferrite.n_components(::MatrixizedInterpolation{vdim1, vdim2}) where {vdim1, vdim2} = vdim1*vdim2
+Ferrite.adjust_dofs_during_distribution(ip::MatrixizedInterpolation) = Ferrite.adjust_dofs_during_distribution(ip.ip)
+
+# Vectorize to reference dimension by default
+function MatrixizedInterpolation(ip::ScalarInterpolation{shape}) where {refdim, shape <: Ferrite.AbstractRefShape{refdim}}
+    return MatrixizedInterpolation{refdim,refdim}(ip)
+end
+
+Base.:(^)(ip::VectorizedInterpolation{vdim1}, vdim2::Int) where {vdim1} = MatrixizedInterpolation{vdim1, vdim2}(ip)
+function Base.literal_pow(::typeof(^), ip::VectorizedInterpolation{vdim1}, ::Val{vdim2}) where {vdim1,vdim2}
+    return MatrixizedInterpolation{vdim1, vdim2}(ip.ip)
+end
+
+function Base.show(io::IO, mime::MIME"text/plain", ip::MatrixizedInterpolation{vdim1, vdim2}) where {vdim1, vdim2}
+    show(io, mime, ip.ip)
+    print(io, "^", vdim1 , "×", vdim2)
+end
+
+# Helper to get number of copies for DoF distribution
+Ferrite.get_n_copies(::MatrixizedInterpolation{vdim1, vdim2}) where {vdim1, vdim2} = vdim1*vdim2
+
+function Ferrite.getnbasefunctions(ipv::MatrixizedInterpolation{vdim1, vdim2}) where {vdim1, vdim2}
+    return vdim1 * vdim2 * getnbasefunctions(ipv.ip)
+end
+function Ferrite.reference_shape_value(ipv::MatrixizedInterpolation{vdim, vdim, shape}, ξ::Tensors.Vec{refdim, T}, I::Int) where {vdim, refdim, shape <: Ferrite.AbstractRefShape{refdim}, T}
+    # First flatten to vector
+    i0, c0 = divrem(I - 1, vdim^2)
+    i = i0 + 1
+    v = Ferrite.reference_shape_value(ipv.ip, ξ, i)
+
+    # Then compute matrix index
+    ci0, cj0 = divrem(c0, vdim)
+    ci = ci0 + 1
+    cj = cj0 + 1
+    return Ferrite.Tensor{2, vdim, T}((k, l) -> k == ci && l == cj ? v : zero(v))
+end
+
+# vdim1 == vdim2 == refdim
+# function shape_gradient_and_value(ipv::MatrixizedInterpolation{dim, dim, shape}, ξ::Vec{dim}, I::Int) where {dim, shape <: AbstractRefShape{dim}}
+#     # TODO order 3 tensor
+#     return invoke(shape_gradient_and_value, Tuple{Interpolation, Vec, Int}, ipv, ξ, I)
+# end
+# vdim1 != vdim2 != refdim
+# function shape_gradient_and_value(ipv::MatrixizedInterpolation{vdim1, vdim2, shape}, ξ::V, I::Int) where {vdim1, vdim2, refdim, shape <: AbstractRefShape{refdim}, T, V <: Vec{refdim, T}}
+#     # Load with dual numbers and compute the value
+#     f = x -> reference_shape_value(ipv, x, I)
+#     ξd = Tensors._load(ξ, Tensors.Tag(f, V))
+#     value_grad = f(ξd)
+#     # Extract the value and gradient
+#     val = Vec{vdim, T}(i -> Tensors.value(value_grad[i]))
+#     grad = zero(MMatrix{vdim, refdim, T})
+#     for (i, vi) in pairs(value_grad)
+#         p = Tensors.partials(vi)
+#         for (j, pj) in pairs(p)
+#             grad[i, j] = pj
+#         end
+#     end
+#     return SArray{Tuple{vdim1, vdim2, }(grad), val
+# end
+
+Ferrite.reference_coordinates(ip::MatrixizedInterpolation) = Ferrite.reference_coordinates(ip.ip)
+
+Ferrite.is_discontinuous(::Type{<:MatrixizedInterpolation{<:Any, <:Any, <:Any, <:Any, ip}}) where {ip} = Ferrite.is_discontinuous(ip)
+
+get_gradient_interpolation(::VectorizedInterpolation{vdim, shape, order, <:Lagrange{shape, order}}) where {sdim,vdim,shape<:Ferrite.AbstractRefShape{sdim},order} = MatrixizedInterpolation{vdim, sdim}(DiscontinuousLagrange{shape, order-1}())
+get_gradient_interpolation_type(::Type{VectorizedInterpolation{vdim, shape, order, <:Lagrange{shape, order}}}) where {sdim,vdim,shape<:Ferrite.AbstractRefShape{sdim},order} = MatrixizedInterpolation{vdim, sdim, shape, order-1, DiscontinuousLagrange{shape, order-1}}
+
+Ferrite.InterpolationInfo(ip::MatrixizedInterpolation) = Ferrite.InterpolationInfo(ip.ip, Ferrite.get_n_copies(ip))
diff --git a/test/runtests.jl b/test/runtests.jl
index e57c23a..1c4e801 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,28 +1,27 @@
 using FerriteViz, Ferrite
 using Test
 
-# _test_tolerance(ip::Interpolation{<:Any,<:Any,1}) = 5e-1
-_test_tolerance(ip::Interpolation) = 1e-6 # Float32 computations are involved!
+# Float32 computations are involved!
+_test_tolerance(ip::Interpolation{<:Any,1}) = 5e-1
+_test_tolerance(ip::Interpolation) = 1e-6
 
-struct MatrixValued <: Ferrite.FieldTrait end
-
-function Ferrite.function_value(::MatrixValued, fe_v::Ferrite.Values{dim}, q_point::Int, u::AbstractVector) where {dim}
+function Ferrite.function_value(fe_v::Ferrite.CellValues{Ferrite.FunctionValues{<:Any, <: FerriteViz.MatrixizedInterpolation{dim,dim}}}, q_point::Int, u::AbstractVector) where {dim}
     n_base_funcs = Ferrite.getn_scalarbasefunctions(fe_v)
     length(u) == n_base_funcs*dim^2 || Ferrite.throw_incompatible_dof_length(length(u), n_base_funcs)
     @boundscheck Ferrite.checkquadpoint(fe_v, q_point)
     val = zero(Tensor{2, dim})
-    
+
     @inbounds for I ∈ 1:n_base_funcs*dim^2
          # First flatten to vector
         i0, c0 = divrem(I - 1, dim^2)
         i = i0 + 1
         v = Ferrite.shape_value(fe_v, q_point, i)
-    
+
         # Then compute matrix index
         ci0, cj0 = divrem(c0, dim)
         ci = ci0 + 1
         cj = cj0 + 1
-        
+
         val += Ferrite.Tensor{2, dim}((k, l) -> k == ci && l == cj ? v*u[I] : zero(v))
     end
 
@@ -32,19 +31,19 @@ end
 @testset "utility operations" begin
     # Check scalar problems
     for (num_elements_per_dim, geo, ip) ∈ [
-                        #    (4,Triangle, Lagrange{2,RefTetrahedron,1}()),
-                           (2,Triangle, Lagrange{2,RefTetrahedron,2}()),
-                           (2,Triangle, Lagrange{2,RefTetrahedron,3}()),
-                        #    (5,Tetrahedron, Lagrange{3,RefTetrahedron,1}()),
-                           (2,Tetrahedron, Lagrange{3,RefTetrahedron,2}()),
-                        #    (4,Quadrilateral, Lagrange{2,RefCube,1}()),
-                           (2,Quadrilateral, Lagrange{2,RefCube,2}()),
-                        #    (4,Hexahedron, Lagrange{3,RefCube,1}()),
-                           (2,Hexahedron, Lagrange{3,RefCube,2}())
+            # (4,Triangle, Lagrange{RefTriangle,1}()),
+            (2,Triangle, Lagrange{RefTriangle,2}()),
+            (2,Triangle, Lagrange{RefTriangle,3}()),
+            # (5,Tetrahedron, Lagrange{RefTetrahedron,1}()),
+            (3,Tetrahedron, Lagrange{RefTetrahedron,2}()),
+            # (4,Quadrilateral, Lagrange{RefQuadrilateral,1}()),
+            (2,Quadrilateral, Lagrange{RefQuadrilateral,2}()),
+            # (4,Hexahedron, Lagrange{RefHexahedron,1}()),
+            (2,Hexahedron, Lagrange{RefHexahedron,2}())
         ]
-        @testset "scalar($num_elements_per_dim, $geo, $ip)" begin
+        @testset failfast=true "scalar($num_elements_per_dim, $geo, $ip)" begin
             # Get solution
-            dim = Ferrite.getdim(ip)
+            dim = Ferrite.getrefdim(ip)
             grid = generate_grid(geo, ntuple(x->num_elements_per_dim, dim));
 
             dh = DofHandler(grid)
@@ -69,10 +68,10 @@ end
 
                 # Check gradient of solution
                 @testset "interpolate_gradient_field" begin
-                    qr = QuadratureRule{dim,Ferrite.getrefshape(ip)}(2) # TODO sample random point
-                    ip_geo = Ferrite.default_interpolation(geo)
+                    qr = QuadratureRule{Ferrite.getrefshape(ip)}(2) # TODO sample random point
+                    ip_geo = Ferrite.geometric_interpolation(geo)
                     ip_grad = Ferrite.getfieldinterpolation(dh_grad, Ferrite.find_field(dh_grad, :gradient))
-                    cellvalues_grad = Ferrite.CellVectorValues(qr, ip_grad, ip_geo)
+                    cellvalues_grad = Ferrite.CellValues(qr, ip_grad, ip_geo)
                     for cell in CellIterator(dh_grad)
                         reinit!(cellvalues_grad, cell)
                         coords = getcoordinates(cell)
@@ -102,13 +101,13 @@ end
             end
         end
 
-        @testset "vector($num_elements_per_dim, $geo, $ip)" begin
+        @testset failfast=true "vector($num_elements_per_dim, $geo, $ip)" begin
             # Get solution
-            dim = Ferrite.getdim(ip)
+            dim = Ferrite.getrefdim(ip)
             grid = generate_grid(geo, ntuple(x->num_elements_per_dim, dim));
 
             dh = DofHandler(grid)
-            add!(dh, :u, dim, ip)
+            add!(dh, :u, ip^dim)
             close!(dh);
 
             # Some test functions with rather complicated gradients
@@ -142,17 +141,17 @@ end
 
                 # Check gradient of solution
                 @testset "interpolate_gradient_field" begin
-                    qr = QuadratureRule{dim,Ferrite.getrefshape(ip)}(2) # TODO sample random point
-                    ip_geo = Ferrite.default_interpolation(geo)
+                    qr = QuadratureRule{Ferrite.getrefshape(ip)}(2) # TODO sample random point
+                    ip_geo = Ferrite.geometric_interpolation(geo)
                     ip_grad = Ferrite.getfieldinterpolation(dh_grad, Ferrite.find_field(dh_grad, :gradient))
-                    cellvalues_grad = Ferrite.CellScalarValues(qr, ip_grad, ip_geo)
+                    cellvalues_grad = Ferrite.CellValues(qr, ip_grad, ip_geo)
                     for cell in CellIterator(dh_grad)
                         reinit!(cellvalues_grad, cell)
                         coords = getcoordinates(cell)
                         uₑ = @views u_grad[celldofs(cell)]
                         for q_point in 1:getnquadpoints(cellvalues_grad)
                             x = spatial_coordinate(cellvalues_grad, q_point, coords)
-                            ∇uₐₚₚᵣₒₓ = function_value(MatrixValued(), cellvalues_grad, q_point, uₑ)
+                            ∇uₐₚₚᵣₒₓ = function_value(cellvalues_grad, q_point, uₑ)
                             ∇uₐₙₐ = Tensors.gradient(f_ana, x)
                             @test all(isapprox.(∇uₐₙₐ, ∇uₐₚₚᵣₒₓ;atol=_test_tolerance(ip)))
                         end
@@ -167,8 +166,7 @@ end
                     for i ∈ 1:size(data_grad, 1)
                         !visible_nodes_grad[i] && continue
                         x = Vec{dim,Float64}(plotter_grad.physical_coords[i])
-                        # Transpose because constructed from Vector and not from Tuple :)
-                        ∇uₐₚₚᵣₒₓ = transpose(Tensor{2,dim,Float64,2*dim}(data_grad[i,:]))
+                        ∇uₐₚₚᵣₒₓ = Tensor{2,dim,Float64,2*dim}(data_grad[i,:])
                         ∇uₐₙₐ = Tensors.gradient(f_ana, x)
                         @test all(isapprox.(∇uₐₙₐ, ∇uₐₚₚᵣₒₓ; atol=_test_tolerance(ip)))
                     end