diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a31ddb90b1..dc5e52ae76 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -13,7 +13,7 @@ jobs: continue-on-error: ${{ matrix.julia-version == 'nightly' }} strategy: matrix: - julia-version: ['1.0', '1', 'nightly'] + julia-version: ['1.6', '^1.7.0-0', 'nightly'] os: ['ubuntu-latest'] include: - os: windows-latest @@ -28,9 +28,7 @@ jobs: - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - if: ${{ matrix.julia-version == '1' && matrix.os == 'ubuntu-latest' }} - uses: codecov/codecov-action@v1 - if: ${{ matrix.julia-version == '1' && matrix.os == 'ubuntu-latest' }} with: file: lcov.info docs: @@ -47,4 +45,5 @@ jobs: env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} + GKSwstype: '100' run: julia --project=docs --color=yes docs/make.jl diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 0000000000..dbcb7ff684 --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,13 @@ +cff-version: 1.2.0 +message: "If you use Ferrite.jl, please cite it as below." +authors: +- family-names: "Carlsson" + given-names: "Kristoffer" +- family-names: "Ekre" + given-names: "Fredrik" +- family-names: "Contributors" + given-names: " " +title: "Ferrite.jl" +version: 0.3.0 +date-released: 2021-03-25 +url: "https://github.com/Ferrite-FEM/Ferrite.jl" diff --git a/Project.toml b/Project.toml index f4bb5e7115..7a8f274c31 100644 --- a/Project.toml +++ b/Project.toml @@ -1,19 +1,21 @@ name = "Ferrite" uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992" -version = "0.3.0" +version = "0.3.1" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] +NearestNeighbors = "0.4" Reexport = "1" Tensors = "1" WriteVTK = "1" -julia = "1" +julia = "1.6" [extras] BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" diff --git a/README.md b/README.md index cad5f62908..1396239efc 100644 --- a/README.md +++ b/README.md @@ -23,5 +23,5 @@ Contributions are very welcome, as are feature requests and suggestions. Please keep in mind that we are part of the Julia community and adhere to the [Julia Community Standards](https://julialang.org/community/standards/). -[docs-dev-img]: https://img.shields.io/badge/docs-dev-blue.svg -[docs-dev-url]: http://ferrite-fem.github.io/Ferrite.jl/dev/ +[docs-dev-img]: https://img.shields.io/badge/docs-latest%20release-blue +[docs-dev-url]: http://ferrite-fem.github.io/Ferrite.jl/ diff --git a/docs/Manifest.toml b/docs/Manifest.toml index ac2b1a17f5..f5a71f2847 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -1,55 +1,84 @@ # This file is machine-generated - editing it directly is not advised +[[ANSIColoredPrinters]] +git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" +uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9" +version = "0.0.1" + [[Adapt]] deps = ["LinearAlgebra"] -git-tree-sha1 = "ffcfa2d345aaee0ef3d8346a073d5dd03c983ebe" +git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.2.0" +version = "3.3.1" + +[[ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" + +[[ArnoldiMethod]] +deps = ["LinearAlgebra", "Random", "StaticArrays"] +git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9" +uuid = "ec485272-7323-5ecc-a04f-4719b315124d" +version = "0.1.0" [[ArrayInterface]] -deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"] -git-tree-sha1 = "ce17bad65d0842b34a15fffc8879a9f68f08a67f" +deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"] +git-tree-sha1 = "b8d49c34c3da35f220e7295659cd0bab8e739fed" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "3.1.6" +version = "3.1.33" [[ArrayLayouts]] deps = ["FillArrays", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "9aa9647b58147a81f7359eacc7d6249ac3a3e3d4" +git-tree-sha1 = "7a92ea1dd16472d18ca1ffcbb7b3cc67d7e78a3f" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" -version = "0.6.4" +version = "0.7.7" [[Artifacts]] -deps = ["Pkg"] -git-tree-sha1 = "c30985d8821e0cd73870b17b0ed0ce6dc44cb744" uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" -version = "1.3.0" [[Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +[[BitTwiddlingConvenienceFunctions]] +deps = ["Static"] +git-tree-sha1 = "652aab0fc0d6d4db4cc726425cadf700e9f473f1" +uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b" +version = "0.1.0" + [[BlockArrays]] deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"] -git-tree-sha1 = "8ecf08b50d71b3fcf3272ee935b3eda653bad5bb" +git-tree-sha1 = "5524e27323cf4c4505699c3fb008c3f772269945" uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" -version = "0.15.2" +version = "0.16.9" [[Bzip2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "c3598e525718abcc440f69cc6d5f60dda0a1b61e" +git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2" uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" -version = "1.0.6+5" +version = "1.0.8+0" + +[[CPUSummary]] +deps = ["Hwloc", "IfElse", "Static"] +git-tree-sha1 = "38d0941d2ce6dd96427fd033d45471e1f26c3865" +uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9" +version = "0.1.5" [[Cairo_jll]] deps = ["Artifacts", "Bzip2_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 = "e2f47f6d8337369411569fd45ae5753ca10394c6" +git-tree-sha1 = "f2202b55d816427cd385a9a4f3ffb226bee80f99" uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a" -version = "1.16.0+6" +version = "1.16.1+0" [[ChainRulesCore]] deps = ["Compat", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "de4f08843c332d355852721adb1592bce7924da3" +git-tree-sha1 = "a325370b9dd0e6bf5656a6f1a7ae80755f8ccc46" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "0.9.29" +version = "1.7.2" + +[[CloseOpenIntervals]] +deps = ["ArrayInterface", "Static"] +git-tree-sha1 = "ce9c0d07ed6e1a4fecd2df6ace144cbd29ba6f37" +uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9" +version = "0.1.2" [[CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] @@ -58,22 +87,27 @@ uuid = "944b1d66-785c-5afd-91f1-9de20f533193" version = "0.7.0" [[ColorSchemes]] -deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"] -git-tree-sha1 = "3141757b5832ee7a0386db87997ee5a23ff20f4d" +deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"] +git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a" uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4" -version = "3.10.2" +version = "3.15.0" [[ColorTypes]] deps = ["FixedPointNumbers", "Random"] -git-tree-sha1 = "32a2b8af383f11cbb65803883837a149d10dfe8a" +git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597" uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" -version = "0.10.12" +version = "0.11.0" [[Colors]] -deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Reexport"] -git-tree-sha1 = "ac5f2213e56ed8a34a3dd2f681f4df1166b34929" +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" -version = "0.12.6" +version = "0.12.8" + +[[CommonSolve]] +git-tree-sha1 = "68a0743f578349ada8bc911a5cbd5a2ef6ed6d1f" +uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" +version = "0.2.0" [[CommonSubexpressions]] deps = ["MacroTools", "Test"] @@ -83,15 +117,19 @@ version = "0.3.0" [[Compat]] deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "919c7f3151e79ff196add81d7f4e45d91bbf420b" +git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "3.25.0" +version = "3.39.0" [[CompilerSupportLibraries_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "8e695f735fca77e9708e795eda62afdb869cbb70" +deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "0.3.4+0" + +[[ConstructionBase]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4" +uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +version = "1.3.0" [[Contour]] deps = ["StaticArrays"] @@ -104,16 +142,22 @@ git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d" uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" version = "4.0.4" +[[DEDataArrays]] +deps = ["ArrayInterface", "DocStringExtensions", "LinearAlgebra", "RecursiveArrayTools", "SciMLBase", "StaticArrays"] +git-tree-sha1 = "31186e61936fbbccb41d809ad4338c9f7addf7ae" +uuid = "754358af-613d-5f8d-9788-280bf1605d4c" +version = "0.2.0" + [[DataAPI]] -git-tree-sha1 = "dfb3b7e89e395be1e25c2ad6d7690dc29cc53b1d" +git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8" uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -version = "1.6.0" +version = "1.9.0" [[DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677" +git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.9" +version = "0.18.10" [[DataValueInterfaces]] git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" @@ -128,6 +172,12 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" deps = ["Mmap"] uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" +[[DiffEqBase]] +deps = ["ArrayInterface", "ChainRulesCore", "DEDataArrays", "DataStructures", "Distributions", "DocStringExtensions", "FastBroadcast", "ForwardDiff", "FunctionWrappers", "IterativeSolvers", "LabelledArrays", "LinearAlgebra", "Logging", "MuladdMacro", "NonlinearSolve", "Parameters", "PreallocationTools", "Printf", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "Requires", "SciMLBase", "Setfield", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "ZygoteRules"] +git-tree-sha1 = "420ad175d5e420e2c55a0ed8a9c18556e6735f80" +uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" +version = "6.73.2" + [[DiffResults]] deps = ["StaticArrays"] git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805" @@ -136,67 +186,105 @@ version = "1.0.3" [[DiffRules]] deps = ["NaNMath", "Random", "SpecialFunctions"] -git-tree-sha1 = "214c3fcac57755cfda163d91c58893a8723f93e9" +git-tree-sha1 = "7220bc21c33e990c14f4a9a319b1d242ebc5b269" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" -version = "1.0.2" +version = "1.3.1" + +[[Distances]] +deps = ["LinearAlgebra", "Statistics", "StatsAPI"] +git-tree-sha1 = "9f46deb4d4ee4494ffb5a40a27a2aced67bdd838" +uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +version = "0.10.4" [[Distributed]] deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" +[[Distributions]] +deps = ["ChainRulesCore", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"] +git-tree-sha1 = "ff7890c74e2eaffbc0b3741811e3816e64b6343d" +uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" +version = "0.25.18" + [[DocStringExtensions]] -deps = ["LibGit2", "Markdown", "Pkg", "Test"] -git-tree-sha1 = "9d4f64f79012636741cf01133158a54b24924c32" +deps = ["LibGit2"] +git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.8.4" +version = "0.8.5" [[Documenter]] -deps = ["Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] -git-tree-sha1 = "3ebb967819b284dc1e3c0422229b58a40a255649" +deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] +git-tree-sha1 = "8b43e37cfb4f4edc2b6180409acc0cebce7fede8" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "0.26.3" +version = "0.27.7" + +[[Downloads]] +deps = ["ArgTools", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" [[EarCut_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "92d8f9f208637e8d2d28c664051a00569c01493d" +git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d" uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5" -version = "2.1.5+1" +version = "2.2.3+0" [[Expat_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "1402e52fcda25064f51c77a9655ce8680b76acf0" +git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f" uuid = "2e619515-83b5-522b-bb60-26c02a35a201" -version = "2.2.7+6" +version = "2.2.10+0" + +[[ExponentialUtilities]] +deps = ["ArrayInterface", "LinearAlgebra", "Printf", "Requires", "SparseArrays"] +git-tree-sha1 = "54b4bd8f88278fd544a566465c943ce4f8da7b7f" +uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18" +version = "1.10.0" + +[[ExprTools]] +git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92" +uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" +version = "0.1.6" [[FFMPEG]] -deps = ["FFMPEG_jll", "x264_jll"] -git-tree-sha1 = "9a73ffdc375be61b0e4516d83d880b265366fe1f" +deps = ["FFMPEG_jll"] +git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8" uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" -version = "0.4.0" +version = "0.4.1" [[FFMPEG_jll]] -deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "LibVPX_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"] -git-tree-sha1 = "3cc57ad0a213808473eafef4845a74766242e05f" +deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"] +git-tree-sha1 = "d8a578692e3077ac998b50c0217dfd67f21d1e5f" uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5" -version = "4.3.1+4" +version = "4.4.0+0" + +[[FastBroadcast]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "26be48918640ce002f5833e8fc537b2ba7ed0234" +uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +version = "0.1.8" + +[[FastClosures]] +git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef" +uuid = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" +version = "0.3.2" [[Ferrite]] -deps = ["LinearAlgebra", "Reexport", "SparseArrays", "Tensors", "WriteVTK"] +deps = ["LinearAlgebra", "NearestNeighbors", "Reexport", "SparseArrays", "Tensors", "WriteVTK"] path = ".." uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992" -version = "0.3.0" +version = "0.3.1" [[FillArrays]] -deps = ["LinearAlgebra", "Random", "SparseArrays"] -git-tree-sha1 = "31939159aeb8ffad1d4d8ee44d07f8558273120a" +deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"] +git-tree-sha1 = "29890dfbc427afa59598b8cfcc10034719bd7744" uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "0.11.7" +version = "0.12.6" [[FiniteDiff]] deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"] -git-tree-sha1 = "f6f80c8f934efd49a286bb5315360be66956dfc4" +git-tree-sha1 = "8b3c09b56acaf3c0e581c66638b85c8650ee9dca" uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" -version = "2.8.0" +version = "2.8.1" [[FixedPointNumbers]] deps = ["Statistics"] @@ -206,9 +294,9 @@ version = "0.8.4" [[Fontconfig_jll]] deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "35895cf184ceaab11fd778b4590144034a167a2f" +git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03" uuid = "a3f928ae-7b40-5064-980b-68af3947d34b" -version = "2.13.1+14" +version = "2.13.93+0" [[Formatting]] deps = ["Printf"] @@ -217,81 +305,125 @@ uuid = "59287772-0a20-5a39-b81b-1366585eb4c0" version = "0.4.2" [[ForwardDiff]] -deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "NaNMath", "Random", "SpecialFunctions", "StaticArrays"] -git-tree-sha1 = "c68fb7481b71519d313114dca639b35262ff105f" +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Printf", "Random", "SpecialFunctions", "StaticArrays"] +git-tree-sha1 = "c4203b60d37059462af370c4f3108fb5d155ff13" uuid = "f6369f11-7733-5829-9624-2563aa707210" -version = "0.10.17" +version = "0.10.20" [[FreeType2_jll]] deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] -git-tree-sha1 = "cbd58c9deb1d304f5a245a0b7eb841a2560cfec6" +git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9" uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7" -version = "2.10.1+5" +version = "2.10.4+0" [[FriBidi_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "0d20aed5b14dd4c9a2453c1b601d08e1149679cc" +git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91" uuid = "559328eb-81f9-559d-9380-de523a88c83c" -version = "1.0.5+6" +version = "1.0.10+0" + +[[FunctionWrappers]] +git-tree-sha1 = "241552bc2209f0fa068b6415b1942cc0aa486bcc" +uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" +version = "1.1.2" + +[[Future]] +deps = ["Random"] +uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" [[GLFW_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"] -git-tree-sha1 = "bd1dbf065d7a4a0bdf7e74dd26cf932dda22b929" +git-tree-sha1 = "dba1e8614e98949abfa60480b13653813d8f0157" uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89" -version = "3.3.3+0" +version = "3.3.5+0" [[GR]] -deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"] -git-tree-sha1 = "12d971c928b7ecf19b748a2c7df6a365690dbf2c" +deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"] +git-tree-sha1 = "c2178cfbc0a5a552e16d097fae508f2024de61a3" uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" -version = "0.55.0" +version = "0.59.0" [[GR_jll]] -deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt_jll", "Zlib_jll", "libpng_jll"] -git-tree-sha1 = "8aee6fa096b0cbdb05e71750c978b96a08c78951" +deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"] +git-tree-sha1 = "cafe0823979a5c9bff86224b3b8de29ea5a44b2e" uuid = "d2c73de3-f751-5644-a686-071e5b155ba9" -version = "0.53.0+0" +version = "0.61.0+0" [[GeometryBasics]] deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"] -git-tree-sha1 = "c7f81b22b6c255861be4007a16bfdeb60e1c7f9b" +git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c" uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326" -version = "0.3.11" +version = "0.4.1" [[Gettext_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"] -git-tree-sha1 = "8c14294a079216000a0bdca5ec5a447f073ddc9d" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"] +git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046" uuid = "78b55507-aeef-58d4-861c-77aaff3498b1" -version = "0.20.1+7" +version = "0.21.0+0" [[Glib_jll]] deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "04690cc5008b38ecbdfede949220bc7d9ba26397" +git-tree-sha1 = "7bf67e9a481712b3dbe9cb3dac852dc4b1162e02" uuid = "7746bdde-850d-59dc-9ae8-88ece973131d" -version = "2.59.0+4" +version = "2.68.3+0" + +[[Graphite2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011" +uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472" +version = "1.3.14+0" [[Grisu]] -git-tree-sha1 = "03d381f65183cb2d0af8b3425fde97263ce9a995" +git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2" uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe" -version = "1.0.0" +version = "1.0.2" [[HTTP]] -deps = ["Base64", "Dates", "IniFile", "MbedTLS", "NetworkOptions", "Sockets", "URIs"] -git-tree-sha1 = "c9f380c76d8aaa1fa7ea9cf97bddbc0d5b15adc2" +deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"] +git-tree-sha1 = "14eece7a3308b4d8be910e265c724a6ba51a9798" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "0.9.5" +version = "0.9.16" + +[[HarfBuzz_jll]] +deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"] +git-tree-sha1 = "8a954fed8ac097d5be04921d595f741115c1b2ad" +uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566" +version = "2.8.1+0" + +[[HostCPUFeatures]] +deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"] +git-tree-sha1 = "3169c8b31863f9a409be1d17693751314241e3eb" +uuid = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0" +version = "0.1.4" + +[[Hwloc]] +deps = ["Hwloc_jll"] +git-tree-sha1 = "92d99146066c5c6888d5a3abc871e6a214388b91" +uuid = "0e44f5e4-bd66-52a0-8798-143a42290a1d" +version = "2.0.0" + +[[Hwloc_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "3395d4d4aeb3c9d31f5929d32760d8baeee88aaf" +uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8" +version = "2.5.0+0" [[IOCapture]] -deps = ["Logging"] -git-tree-sha1 = "377252859f740c217b936cebcd918a44f9b53b59" +deps = ["Logging", "Random"] +git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a" uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" -version = "0.1.1" +version = "0.2.2" [[IfElse]] git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef" uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" version = "0.1.0" +[[Inflate]] +git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c" +uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" +version = "0.1.2" + [[IniFile]] deps = ["Test"] git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8" @@ -302,6 +434,11 @@ version = "0.5.0" deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +[[IrrationalConstants]] +git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94" +uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" +version = "0.1.0" + [[IterTools]] git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18" uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e" @@ -309,9 +446,9 @@ version = "1.3.0" [[IterativeSolvers]] deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"] -git-tree-sha1 = "6f5ef3206d9dc6510a8b8e2334b96454a2ade590" +git-tree-sha1 = "1a8c6237e78b714e901e406c096fc8a65528af7d" uuid = "42fd0dbc-a981-5370-80f2-aaf504508153" -version = "0.9.0" +version = "0.9.1" [[IteratorInterfaceExtensions]] git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856" @@ -319,21 +456,22 @@ uuid = "82899510-4779-5014-852e-03e436cf321d" version = "1.0.0" [[JLLWrappers]] -git-tree-sha1 = "a431f5f2ca3f4feef3bd7a5e94b8b8d4f2f647a0" +deps = ["Preferences"] +git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e" uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.2.0" +version = "1.3.0" [[JSON]] deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4" +git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.1" +version = "0.21.2" [[JpegTurbo_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "9aff0587d9603ea0de2c6f6300d9f9492bbefbd3" +git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943" uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8" -version = "2.0.1+3" +version = "2.1.0+0" [[KrylovMethods]] deps = ["LinearAlgebra", "Printf", "SparseArrays"] @@ -343,51 +481,69 @@ version = "0.6.0" [[LAME_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "df381151e871f41ee86cee4f5f6fd598b8a68826" +git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c" uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d" -version = "3.100.0+3" +version = "3.100.1+0" [[LZO_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "f128cd6cd05ffd6d3df0523ed99b90ff6f9b349a" +git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6" uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac" -version = "2.10.0+3" +version = "2.10.1+0" [[LaTeXStrings]] git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104" uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" version = "1.2.1" +[[LabelledArrays]] +deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"] +git-tree-sha1 = "8f5fd068dfee92655b79e0859ecad8b492dfe8b1" +uuid = "2ee39098-c373-598a-b85f-a56591580800" +version = "1.6.5" + [[Latexify]] deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"] -git-tree-sha1 = "fbc08b5a78e264ba3d19da90b36ce1789ca67a40" +git-tree-sha1 = "a4b12a1bd2ebade87891ab7e36fdbce582301a92" uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" -version = "0.14.11" +version = "0.15.6" + +[[LayoutPointers]] +deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static"] +git-tree-sha1 = "d2bda6aa0b03ce6f141a2dc73d0bcb7070131adc" +uuid = "10f19ff3-798f-405d-979b-55457f8fc047" +version = "0.1.3" + +[[LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" + +[[LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" [[LibGit2]] -deps = ["Printf"] +deps = ["Base64", "NetworkOptions", "Printf", "SHA"] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" -[[LibVPX_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "85fcc80c3052be96619affa2fe2e6d2da3908e11" -uuid = "dd192d2f-8180-539f-9fb4-cc70b1dcf69a" -version = "1.9.0+1" +[[LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" [[Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" [[Libffi_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "a2cd088a88c0d37eef7d209fd3d8712febce0d90" +git-tree-sha1 = "761a393aeccd6aa92ec3515e428c26bf99575b3b" uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490" -version = "3.2.1+4" +version = "3.2.2+0" [[Libgcrypt_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"] -git-tree-sha1 = "b391a18ab1170a2e568f9fb8d83bc7c780cb9999" +git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae" uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4" -version = "1.8.5+4" +version = "1.8.7+0" [[Libglvnd_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"] @@ -397,33 +553,39 @@ version = "1.3.0+3" [[Libgpg_error_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "ec7f2e8ad5c9fa99fc773376cdbc86d9a5a23cb7" +git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9" uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8" -version = "1.36.0+3" +version = "1.42.0+0" [[Libiconv_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "8e924324b2e9275a51407a4e06deb3455b1e359f" +git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778" uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" -version = "1.16.0+7" +version = "1.16.1+1" [[Libmount_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "51ad0c01c94c1ce48d5cad629425035ad030bfd5" +git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73" uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9" -version = "2.34.0+3" +version = "2.35.0+0" [[Libtiff_jll]] deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"] -git-tree-sha1 = "291dd857901f94d683973cdf679984cdf73b56d0" +git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9" uuid = "89763e89-9b03-5906-acba-b20f662cd828" -version = "4.1.0+2" +version = "4.3.0+0" [[Libuuid_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "f879ae9edbaa2c74c922e8b85bb83cc84ea1450b" +git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066" uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" -version = "2.34.0+7" +version = "2.36.0+0" + +[[LightGraphs]] +deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] +git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6" +uuid = "093fc24a-ae57-5d10-9952-331d41423f4d" +version = "1.3.5" [[LightXML]] deps = ["Libdl", "XML2_jll"] @@ -442,19 +604,36 @@ deps = ["Libdl"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[Literate]] -deps = ["Base64", "JSON", "REPL"] -git-tree-sha1 = "32b517d4d8219d3bbab199de3416ace45010bdb3" +deps = ["Base64", "IOCapture", "JSON", "REPL"] +git-tree-sha1 = "bbebc3c14dbfbe76bfcbabf0937481ac84dc86ef" uuid = "98b081ad-f1c9-55d3-8b20-4c87d4299306" -version = "2.8.0" +version = "2.9.3" + +[[LogExpFunctions]] +deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "34dc30f868e368f8a17b728a1238f3fcda43931a" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "0.3.3" [[Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" +[[LoopVectorization]] +deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "Requires", "SIMDDualNumbers", "SLEEFPirates", "Static", "ThreadingUtilities", "UnPack", "VectorizationBase"] +git-tree-sha1 = "b97b8bcf653acaae65a9b57a3a129786c98a0907" +uuid = "bdcacae8-1622-11e9-2a5c-532679323890" +version = "0.12.84" + [[MacroTools]] deps = ["Markdown", "Random"] -git-tree-sha1 = "6a8a2a625ab0dea913aba95c11370589e0239ff0" +git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.6" +version = "0.5.8" + +[[ManualMemory]] +git-tree-sha1 = "9cb207b18148b2199db259adfa923b45593fe08e" +uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667" +version = "0.1.6" [[Markdown]] deps = ["Base64"] @@ -467,10 +646,8 @@ uuid = "739be429-bea8-5141-9913-cc70e7f3736d" version = "1.0.3" [[MbedTLS_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "0eef589dd1c26a3ac9d753fe1a8bcad63f956fa6" +deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.16.8+1" [[Measures]] git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f" @@ -479,90 +656,136 @@ version = "0.3.1" [[Missings]] deps = ["DataAPI"] -git-tree-sha1 = "f8c673ccc215eb50fcadb285f522420e29e69e1c" +git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f" uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" -version = "0.4.5" +version = "1.0.2" [[Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" +[[MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" + +[[MuladdMacro]] +git-tree-sha1 = "c6190f9a7fc5d9d5915ab29f2134421b12d24a68" +uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +version = "0.2.2" + [[NLSolversBase]] deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] -git-tree-sha1 = "50608f411a1e178e0129eab4110bd56efd08816f" +git-tree-sha1 = "144bab5b1443545bc4e791536c9f1eacb4eed06a" uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" -version = "7.8.0" +version = "7.8.1" + +[[NLsolve]] +deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"] +git-tree-sha1 = "019f12e9a1a7880459d0173c182e6a99365d7ac1" +uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +version = "4.5.1" [[NaNMath]] git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" version = "0.3.5" +[[NearestNeighbors]] +deps = ["Distances", "StaticArrays"] +git-tree-sha1 = "16baacfdc8758bc374882566c9187e785e85c2f0" +uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce" +version = "0.4.9" + [[NetworkOptions]] -git-tree-sha1 = "ed3157f48a05543cce9b241e1f2815f7e843d96e" uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" -version = "1.2.0" + +[[NonlinearSolve]] +deps = ["ArrayInterface", "FiniteDiff", "ForwardDiff", "IterativeSolvers", "LinearAlgebra", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "SciMLBase", "Setfield", "StaticArrays", "UnPack"] +git-tree-sha1 = "e9ffc92217b8709e0cf7b8808f6223a4a0936c95" +uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +version = "0.3.11" + +[[OffsetArrays]] +deps = ["Adapt"] +git-tree-sha1 = "c0e9e582987d36d5a61e650e6e543b9e44d9914b" +uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +version = "1.10.7" [[Ogg_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "a42c0f138b9ebe8b58eba2271c5053773bde52d0" +git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f" uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051" -version = "1.3.4+2" +version = "1.3.5+0" + +[[OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" [[OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "71bbbc616a1d710879f5a1021bcba65ffba6ce58" +git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "1.1.1+6" +version = "1.1.10+0" [[OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "9db77584158d0ab52307f8c04f8e7c08ca76b5b3" +git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" -version = "0.5.3+4" +version = "0.5.5+0" [[Optim]] deps = ["Compat", "FillArrays", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"] -git-tree-sha1 = "d34366a3abc25c41f88820762ef7dfdfe9306711" +git-tree-sha1 = "7863df65dbb2a0fa8f85fcaf0a41167640d2ebed" uuid = "429524aa-4258-5aef-a3af-852621145aeb" -version = "1.3.0" +version = "1.4.1" [[Opus_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "f9d57f4126c39565e05a2b0264df99f497fc6f37" +git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720" uuid = "91d4177d-7536-5919-b921-800302f37372" -version = "1.3.1+3" +version = "1.3.2+0" [[OrderedCollections]] -git-tree-sha1 = "4fa2ba51070ec13fcc7517db714445b4ab986bdf" +git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.4.0" +version = "1.4.1" + +[[OrdinaryDiffEq]] +deps = ["Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "ExponentialUtilities", "FastClosures", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "Logging", "LoopVectorization", "MacroTools", "MuladdMacro", "NLsolve", "Polyester", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"] +git-tree-sha1 = "4341419e2badc4efd259bfd67e0726329c454ef0" +uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +version = "5.64.1" [[PCRE_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "1b556ad51dceefdbf30e86ffa8f528b73c7df2bb" +git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488" uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc" -version = "8.42.0+4" +version = "8.44.0+0" + +[[PDMats]] +deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] +git-tree-sha1 = "4dd403333bcf0909341cfe57ec115152f937d7d8" +uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" +version = "0.11.1" [[Parameters]] deps = ["OrderedCollections", "UnPack"] -git-tree-sha1 = "2276ac65f1e236e0a6ea70baff3f62ad4c625345" +git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe" uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a" -version = "0.12.2" +version = "0.12.3" [[Parsers]] deps = ["Dates"] -git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc" +git-tree-sha1 = "a8709b968a1ea6abc2dc1967cb1db6ac9a00dfb6" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "1.1.0" +version = "2.0.5" [[Pixman_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "6a20a83c1ae86416f0a5de605eaea08a552844a3" +git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29" uuid = "30392449-352a-5448-841d-b1acce4e97dc" -version = "0.40.0+0" +version = "0.40.1+0" [[Pkg]] -deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] +deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [[PlotThemes]] @@ -573,15 +796,27 @@ version = "2.0.1" [[PlotUtils]] deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"] -git-tree-sha1 = "ae9a295ac761f64d8c2ec7f9f24d21eb4ffba34d" +git-tree-sha1 = "b084324b4af5a438cd63619fd006614b3b20b87b" uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043" -version = "1.0.10" +version = "1.0.15" [[Plots]] -deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"] -git-tree-sha1 = "40031283dbb5ae602aa132f423daedfc18c82069" +deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"] +git-tree-sha1 = "6841db754bd01a91d281370d9a0f8787e220ae08" uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -version = "1.11.0" +version = "1.22.4" + +[[Polyester]] +deps = ["ArrayInterface", "BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "ManualMemory", "PolyesterWeave", "Requires", "Static", "StrideArraysCore", "ThreadingUtilities"] +git-tree-sha1 = "97794179584fbb0336821d6c03c93682f19803bf" +uuid = "f517fe37-dbe3-4b94-8317-1923a5111588" +version = "0.5.3" + +[[PolyesterWeave]] +deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"] +git-tree-sha1 = "371a19bb801c1b420b29141750f3a34d6c6634b9" +uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad" +version = "0.1.0" [[PositiveFactorizations]] deps = ["LinearAlgebra"] @@ -589,24 +824,42 @@ git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20" uuid = "85a6dd25-e78a-55b7-8502-1745935b8125" version = "0.2.4" +[[PreallocationTools]] +deps = ["ArrayInterface", "ForwardDiff", "LabelledArrays"] +git-tree-sha1 = "361c1f60ffdeeddf02f36b463ab8b138194e5f25" +uuid = "d236fae5-4411-538c-8e31-a6e3d9e00b46" +version = "0.1.1" + +[[Preferences]] +deps = ["TOML"] +git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.2.2" + [[Printf]] deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" [[ProgressMeter]] deps = ["Distributed", "Printf"] -git-tree-sha1 = "6e9c89cba09f6ef134b00e10625590746ba1e036" +git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8" uuid = "92933f4c-e287-5a05-a399-4b506db050ca" -version = "1.5.0" +version = "1.7.1" -[[Qt_jll]] +[[Qt5Base_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"] -git-tree-sha1 = "588b799506f0b956a7e6e18f767dfa03cf333f26" -uuid = "ede63266-ebff-546c-83e0-1c6fb6d0efc8" -version = "5.15.2+2" +git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8" +uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1" +version = "5.15.3+0" + +[[QuadGK]] +deps = ["DataStructures", "LinearAlgebra"] +git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39" +uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +version = "2.4.2" [[REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets"] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" [[Random]] @@ -614,20 +867,32 @@ deps = ["Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [[RecipesBase]] -git-tree-sha1 = "b3fb709f3c97bfc6e948be68beeecb55a0b340ae" +git-tree-sha1 = "44a75aa7a527910ee3d1751d1f0e4148698add9e" uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" -version = "1.1.1" +version = "1.1.2" [[RecipesPipeline]] deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"] -git-tree-sha1 = "c4d54a78e287de7ec73bbc928ce5eb3c60f80b24" +git-tree-sha1 = "7ad0dfa8d03b7bcf8c597f59f5292801730c55b8" uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c" -version = "0.3.1" +version = "0.4.1" + +[[RecursiveArrayTools]] +deps = ["ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"] +git-tree-sha1 = "ff7495c78a192ff7d59531d9f14db300c847a4bc" +uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" +version = "2.19.1" + +[[RecursiveFactorization]] +deps = ["LinearAlgebra", "LoopVectorization", "Polyester", "StrideArraysCore", "TriangularSolve"] +git-tree-sha1 = "575c18c6b00ce409f75d96fefe33ebe01575457a" +uuid = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" +version = "0.2.4" [[Reexport]] -git-tree-sha1 = "57d8440b0c7d98fc4f889e478e80f268d534c9d5" +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" uuid = "189a3867-3050-52da-a836-e630ba90ab69" -version = "1.0.0" +version = "1.2.2" [[Requires]] deps = ["UUIDs"] @@ -635,118 +900,218 @@ git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621" uuid = "ae029012-a4dd-5104-9daa-d747884805df" version = "1.1.3" +[[Rmath]] +deps = ["Random", "Rmath_jll"] +git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f" +uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa" +version = "0.7.0" + +[[Rmath_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7" +uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f" +version = "0.3.0+0" + [[SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" [[SIMD]] -git-tree-sha1 = "b2202f63f62d8c1214ed16c18319f72fa74adc4d" +git-tree-sha1 = "9ba33637b24341aba594a2783a502760aa0bff04" uuid = "fdea26ae-647d-5447-a871-4b548cad5224" -version = "3.2.1" +version = "3.3.1" + +[[SIMDDualNumbers]] +deps = ["ForwardDiff", "IfElse", "SLEEFPirates", "VectorizationBase"] +git-tree-sha1 = "62c2da6eb66de8bb88081d20528647140d4daa0e" +uuid = "3cdde19b-5bb0-4aaf-8931-af3e248e098b" +version = "0.1.0" + +[[SIMDTypes]] +git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c" +uuid = "94e857df-77ce-4151-89e5-788b33177be4" +version = "0.1.0" + +[[SLEEFPirates]] +deps = ["IfElse", "Static", "VectorizationBase"] +git-tree-sha1 = "2e8150c7d2a14ac68537c7aac25faa6577aff046" +uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa" +version = "0.6.27" + +[[SciMLBase]] +deps = ["ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "RecipesBase", "RecursiveArrayTools", "StaticArrays", "Statistics", "Tables", "TreeViews"] +git-tree-sha1 = "2577b69c151374cc505543b1ee3df98b1ac3abb0" +uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +version = "1.19.1" [[Scratch]] deps = ["Dates"] -git-tree-sha1 = "ad4b278adb62d185bbcb6864dc24959ab0627bf6" +git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda" uuid = "6c6a2e73-6563-6170-7368-637461726353" -version = "1.0.3" +version = "1.1.0" [[Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" +[[Setfield]] +deps = ["ConstructionBase", "Future", "MacroTools", "Requires"] +git-tree-sha1 = "fca29e68c5062722b5b4435594c3d1ba557072a3" +uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" +version = "0.7.1" + [[SharedArrays]] deps = ["Distributed", "Mmap", "Random", "Serialization"] uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" [[Showoff]] deps = ["Dates", "Grisu"] -git-tree-sha1 = "ee010d8f103468309b8afac4abb9be2e18ff1182" +git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de" uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f" -version = "0.3.2" +version = "1.0.3" + +[[SimpleTraits]] +deps = ["InteractiveUtils", "MacroTools"] +git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231" +uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d" +version = "0.9.4" [[Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" [[SortingAlgorithms]] -deps = ["DataStructures", "Random", "Test"] -git-tree-sha1 = "03f5898c9959f8115e30bc7226ada7d0df554ddd" +deps = ["DataStructures"] +git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508" uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" -version = "0.3.1" +version = "1.0.1" [[SparseArrays]] deps = ["LinearAlgebra", "Random"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +[[SparseDiffTools]] +deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "LightGraphs", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays", "VertexSafeGraphs"] +git-tree-sha1 = "36a4d27a02af48a1eafd2baff58b32deeeb68926" +uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804" +version = "1.16.5" + [[SpecialFunctions]] -deps = ["ChainRulesCore", "OpenSpecFun_jll"] -git-tree-sha1 = "5919936c0e92cff40e57d0ddf0ceb667d42e5902" +deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "793793f1df98e3d7d554b65a107e9c9a6399a6ed" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "1.3.0" +version = "1.7.0" [[Static]] deps = ["IfElse"] -git-tree-sha1 = "ddec5466a1d2d7e58adf9a427ba69763661aacf6" +git-tree-sha1 = "a8f30abc7c64a39d389680b74e749cf33f872a70" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.2.4" +version = "0.3.3" [[StaticArrays]] deps = ["LinearAlgebra", "Random", "Statistics"] -git-tree-sha1 = "9da72ed50e94dbff92036da395275ed114e04d49" +git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.0.1" +version = "1.2.13" [[Statistics]] deps = ["LinearAlgebra", "SparseArrays"] uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +[[StatsAPI]] +git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510" +uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" +version = "1.0.0" + [[StatsBase]] -deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"] -git-tree-sha1 = "a83fa3021ac4c5a918582ec4721bc0cf70b495a9" +deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] +git-tree-sha1 = "65fb73045d0e9aaa39ea9a29a5e7506d9ef6511f" uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.33.4" +version = "0.33.11" + +[[StatsFuns]] +deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] +git-tree-sha1 = "95072ef1a22b057b1e80f73c2a89ad238ae4cfff" +uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" +version = "0.9.12" + +[[StrideArraysCore]] +deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "ManualMemory", "Requires", "SIMDTypes", "Static", "ThreadingUtilities"] +git-tree-sha1 = "1258e25e171aec339866f283a11e7d75867e77d7" +uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da" +version = "0.2.4" [[StructArrays]] -deps = ["Adapt", "DataAPI", "Tables"] -git-tree-sha1 = "26ea43b4be7e919a2390c3c0f824e7eb4fc19a0a" +deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"] +git-tree-sha1 = "2ce41e0d042c60ecd131e9fb7154a3bfadbf50d3" uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" -version = "0.5.0" +version = "0.6.3" + +[[SuiteSparse]] +deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] +uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" + +[[TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" [[TableTraits]] deps = ["IteratorInterfaceExtensions"] -git-tree-sha1 = "b1ad568ba658d8cbb3b892ed5380a6f3e781a81e" +git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39" uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" -version = "1.0.0" +version = "1.0.1" [[Tables]] deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"] -git-tree-sha1 = "a9ff3dfec713c6677af435d6a6d65f9744feef67" +git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -version = "1.4.1" +version = "1.6.0" + +[[Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" [[Tensors]] deps = ["ForwardDiff", "LinearAlgebra", "SIMD", "Statistics"] -git-tree-sha1 = "778b3ee366af2cd6082e2c072b7e3766e2735f8e" +git-tree-sha1 = "46662e28c5d32a77f21a57d2263aa3bf8d1edee4" uuid = "48a634ad-e948-5137-8d70-aa71f2a747f4" -version = "1.4.4" +version = "1.6.1" [[Test]] -deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +[[ThreadingUtilities]] +deps = ["ManualMemory"] +git-tree-sha1 = "03013c6ae7f1824131b2ae2fc1d49793b51e8394" +uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5" +version = "0.4.6" + [[TimerOutputs]] -deps = ["Printf"] -git-tree-sha1 = "32cdbe6cd2d214c25a0b88f985c9e0092877c236" +deps = ["ExprTools", "Printf"] +git-tree-sha1 = "7cb456f358e8f9d102a8b25e8dfedf58fa5689bc" uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" -version = "0.5.8" +version = "0.5.13" [[TranscodingStreams]] deps = ["Random", "Test"] -git-tree-sha1 = "7c53c35547de1c5b9d46a4797cf6d8253807108c" +git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" -version = "0.9.5" +version = "0.9.6" + +[[TreeViews]] +deps = ["Test"] +git-tree-sha1 = "8d0d7a3fe2f30d6a7f833a5f19f7c7a5b396eae6" +uuid = "a2a6695c-b41b-5b7d-aed9-dbfdeacea5d7" +version = "0.3.0" + +[[TriangularSolve]] +deps = ["CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "LoopVectorization", "Polyester", "Static", "VectorizationBase"] +git-tree-sha1 = "ed55426a514db35f58d36c3812aae89cfc057401" +uuid = "d5829a12-d9aa-46ab-831f-fb7c9ab06edf" +version = "0.1.6" [[URIs]] -git-tree-sha1 = "7855809b88d7b16e9b029afd17880930626f54a2" +git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355" uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4" -version = "1.2.0" +version = "1.3.0" [[UUIDs]] deps = ["Random", "SHA"] @@ -762,15 +1127,27 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" [[UnicodePlots]] deps = ["Crayons", "Dates", "SparseArrays", "StatsBase"] -git-tree-sha1 = "1a63e6eea76b291378ff9f95801f8b6d96213208" +git-tree-sha1 = "f1d09f14722f5f3cef029bcb031be91a92613ae9" uuid = "b8865327-cd53-5732-bb35-84acbb429228" -version = "1.3.0" +version = "2.4.6" + +[[VectorizationBase]] +deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "Hwloc", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static"] +git-tree-sha1 = "a6ca373834ec22aa850c09392399add8f13d476a" +uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f" +version = "0.21.13" + +[[VertexSafeGraphs]] +deps = ["LightGraphs"] +git-tree-sha1 = "b9b450c99a3ca1cc1c6836f560d8d887bcbe356e" +uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" +version = "0.1.2" [[Wayland_jll]] deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"] -git-tree-sha1 = "dc643a9b774da1c2781413fd7b6dcd2c56bb8056" +git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23" uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89" -version = "1.17.0+4" +version = "1.19.0+0" [[Wayland_protocols_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll"] @@ -779,22 +1156,22 @@ uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91" version = "1.18.0+4" [[WriteVTK]] -deps = ["Base64", "CodecZlib", "FillArrays", "LightXML", "Random", "TranscodingStreams"] -git-tree-sha1 = "37eef911a9c5211e0ae4362dc0477cfe6c537ffa" +deps = ["Base64", "CodecZlib", "FillArrays", "LightXML", "TranscodingStreams"] +git-tree-sha1 = "c3403143cecb391ea51fc133be82b024e4ce720b" uuid = "64499a7a-5c06-52f2-abe2-ccb03c286192" -version = "1.9.1" +version = "1.11.0" [[XML2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "be0db24f70aae7e2b89f2f3092e93b8606d659a6" +git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a" uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" -version = "2.9.10+3" +version = "2.9.12+0" [[XSLT_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Pkg", "XML2_jll"] -git-tree-sha1 = "2b3eac39df218762d2d005702d601cd44c997497" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"] +git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a" uuid = "aed1982a-8fda-507f-9586-7b0439959a61" -version = "1.1.33+4" +version = "1.1.34+0" [[Xorg_libX11_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"] @@ -923,52 +1300,64 @@ uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10" version = "1.4.0+3" [[Zlib_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "320228915c8debb12cb434c59057290f0834dbf6" +deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.11+18" [[Zstd_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "2c1332c54931e83f8f94d310fa447fd743e8d600" +git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6" uuid = "3161d3a3-bdf6-5164-811a-617609db77b4" -version = "1.4.8+0" +version = "1.5.0+0" + +[[ZygoteRules]] +deps = ["MacroTools"] +git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0" +uuid = "700de1a5-db45-46bc-99cf-38207098b444" +version = "0.2.2" [[libass_jll]] -deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] -git-tree-sha1 = "acc685bcf777b2202a904cdcb49ad34c2fa1880c" +deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] +git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47" uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0" -version = "0.14.0+4" +version = "0.15.1+0" [[libfdk_aac_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "7a5780a0d9c6864184b3a2eeeb833a0c871f00ab" +git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55" uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280" -version = "0.1.6+4" +version = "2.0.2+0" [[libpng_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] -git-tree-sha1 = "6abbc424248097d69c0c87ba50fcb0753f93e0ee" +git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c" uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f" -version = "1.6.37+6" +version = "1.6.38+0" [[libvorbis_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"] -git-tree-sha1 = "fa14ac25af7a4b8a7f61b287a124df7aab601bcd" +git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb" uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a" -version = "1.3.6+6" +version = "1.3.7+0" + +[[nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" + +[[p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" [[x264_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "d713c1ce4deac133e3334ee12f4adff07f81778f" +git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2" uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a" -version = "2020.7.14+2" +version = "2021.5.5+0" [[x265_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "487da2f8f2f0c8ee0e83f39d13037d6bbf0a45ab" +git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9" uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76" -version = "3.0.0+3" +version = "3.5.0+0" [[xkbcommon_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"] diff --git a/docs/Project.toml b/docs/Project.toml index 85e0b76e10..e970d2d1d6 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,9 +8,11 @@ KrylovMethods = "9a2cd570-f05c-5dc1-9209-93ad6f5727f7" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Optim = "429524aa-4258-5aef-a3af-852621145aeb" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" diff --git a/docs/generate.jl b/docs/generate.jl index 836dfd2ef2..cc8aa37737 100644 --- a/docs/generate.jl +++ b/docs/generate.jl @@ -9,7 +9,12 @@ for example in readdir(EXAMPLEDIR) input = abspath(joinpath(EXAMPLEDIR, example)) script = Literate.script(input, GENERATEDDIR) code = strip(read(script, String)) - mdpost(str) = replace(str, "@__CODE__" => code) + + # remove "hidden" lines which are not shown in the markdown + line_ending_symbol = occursin(code, "\r\n") ? "\r\n" : "\n" + code_clean = join(filter(x->!endswith(x,"#hide"),split(code, r"\n|\r\n")), line_ending_symbol) + + mdpost(str) = replace(str, "@__CODE__" => code_clean) Literate.markdown(input, GENERATEDDIR, postprocess = mdpost) Literate.notebook(input, GENERATEDDIR, execute = true) elseif any(endswith.(example, [".png", ".jpg", ".gif"])) diff --git a/docs/make.jl b/docs/make.jl index f3af6a8f6e..eed33f52eb 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,7 +7,7 @@ include("generate.jl") GENERATEDEXAMPLES = [joinpath("examples", f) for f in ( "heat_equation.md", - "l2_projection.md", + "postprocessing.md", "helmholtz.md", "incompressible_elasticity.md", "hyperelasticity.md", @@ -16,7 +16,8 @@ GENERATEDEXAMPLES = [joinpath("examples", f) for f in ( "transient_heat_equation.md", "landau.md", "linear_shell.md", - "quasi_incompressible_hyperelasticity.md" + "quasi_incompressible_hyperelasticity.md", + "ns_vs_diffeq.md", )] # Build documentation. diff --git a/docs/src/index.md b/docs/src/index.md index cc6d670e96..6a5e598cb7 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -31,13 +31,13 @@ You can install Ferrite from the Pkg REPL (press `]` in the Julia REPL to enter `pkg>` mode): ``` -pkg> add https://github.com/Ferrite-FEM/Ferrite.jl.git +pkg> add Ferrite ``` !!! note Alternative installation method: ```julia - julia> import Pkg; Pkg.add(PackageSpec(url = "https://github.com/Ferrite-FEM/Ferrite.jl.git")) + julia> import Pkg; Pkg.add("Ferrite") ``` To load the package, use diff --git a/docs/src/literate/heat_square_pointevaluation.png b/docs/src/literate/heat_square_pointevaluation.png new file mode 100644 index 0000000000..2d7763ddef Binary files /dev/null and b/docs/src/literate/heat_square_pointevaluation.png differ diff --git a/docs/src/literate/l2_projection.jl b/docs/src/literate/l2_projection.jl deleted file mode 100644 index 2acb081389..0000000000 --- a/docs/src/literate/l2_projection.jl +++ /dev/null @@ -1,85 +0,0 @@ -# # L2-projection -# -# ![](heat_square_fluxes.png) - - -# ## Introduction -# -# This example continues from the Heat equation example, where the temperature field was -# determined on a square domain. In this example, we first compute the heat flux in each -# integration point (based on the solved temperature field) and then we do an L2-projection -# of the fluxes to the nodes of the mesh. By doing this, we can more easily visualize -# integration points quantities. -# -# The L2-projection is defined as follows: Find projection ``q(\boldsymbol{x}) \in L_2(\Omega)`` such that -# ```math -# \int v q \ \mathrm{d}\Omega = \int v d \ \mathrm{d}\Omega \quad \forall v \in L_2(\Omega), -# ``` -# where ``d`` is the quadrature data to project. Since the flux is a vector the projection function -# will be solved with multiple right hand sides, e.g. with ``d = q_x`` and ``d = q_y`` for this 2D problem. -# -# Ferrite has functionality for doing much of this automatically, as displayed in the code below. -# In particular [`L2Projector`](@ref) for assembling the left hand side, and -# [`project`](@ref) for assembling the right hand sides and solving for the projection. - -# ## Implementation -# -# Start by simply running the Heat equation example to solve the problem -include("heat_equation.jl"); - - -# Next we define a function that computes the heat flux for each integration point in the domain. -# Fourier's law is adopted, where the conductivity tensor is assumed to be isotropic with unit -# conductivity ``\lambda = 1 ⇒ q = - \nabla u``, where ``u`` is the temperature. -function compute_heat_fluxes(cellvalues::CellScalarValues{dim,T}, dh::DofHandler, a) where {dim,T} - - n = getnbasefunctions(cellvalues) - cell_dofs = zeros(Int, n) - nqp = getnquadpoints(cellvalues) - - ## Allocate storage for the fluxes to store - q = [Vec{2,T}[] for _ in 1:getncells(dh.grid)] - - for (cell_num, cell) in enumerate(CellIterator(dh)) - q_cell = q[cell_num] - celldofs!(cell_dofs, dh, cell_num) - aᵉ = a[cell_dofs] - reinit!(cellvalues, cell) - - for q_point in 1:nqp - q_qp = - function_gradient(cellvalues, q_point, aᵉ) - push!(q_cell, q_qp) - end - end - return q -end -#md nothing # hide - -# Now call the function to get all the fluxes. -q_gp = compute_heat_fluxes(cellvalues, dh, u); - -# Next, create an [`L2Projector`](@ref) using the same interpolation as was used to approximate the -# temperature field. On instantiation, the projector assembles the coefficient matrix `M` and -# computes the Cholesky factorization of it. By doing so, the projector can be reused without -# having to invert `M` every time. -projector = L2Projector(ip, grid); - -# Project the integration point values to the nodal values -q_nodes = project(projector, q_gp, qr); - - -# ## Exporting to VTK -# To visualize the heat flux, we export the projected field `q_nodes` -# to a VTK-file, which can be viewed in e.g. [ParaView](https://www.paraview.org/). -vtk_grid("heat_equation_flux", grid) do vtk - vtk_point_data(vtk, q_nodes, "q") -end; - -#md # ## [Plain Program](@id l2_projection-plain-program) -#md # -#md # Below follows a version of the program without any comments. -#md # The file is also available here: [l2_projection.jl](l2_projection.jl) -#md # -#md # ```julia -#md # @__CODE__ -#md # ``` diff --git a/docs/src/literate/ns_vs_diffeq.jl b/docs/src/literate/ns_vs_diffeq.jl new file mode 100644 index 0000000000..2317173d94 --- /dev/null +++ b/docs/src/literate/ns_vs_diffeq.jl @@ -0,0 +1,500 @@ +# # Incompressible Navier-Stokes Equations via DifferentialEquations.jl +# +# ![](https://user-images.githubusercontent.com/9196588/134514213-76d91d34-19ab-47c2-957e-16bb0c8669e1.gif) +# +# +# In this example we focus on a simple but visually appealing problem from +# fluid dynamics, namely vortex shedding. This problem is also known as +# [von-Karman vortex streets](https://en.wikipedia.org/wiki/K%C3%A1rm%C3%A1n_vortex_street). Within this example, we show how to utilize [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl) +# in tandem with Ferrite.jl to solve this space-time problem. To keep things simple we use a naive approach +# to discretize the system. +# +# ## Remarks on DifferentialEquations.jl +# +# Many "time step solvers" of [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl) assume that that the +# problem is provided in mass matrix form. The incompressible Navier-Stokes +# equations as stated above yield a DAE in this form after applying a spatial +# discretization technique - in our case FEM. The mass matrix form of ODEs and DAEs +# is given as: +# ```math +# M(t) \mathrm{d}_t u = f(u,t) +# ``` +# where $M$ is a possibly time-dependent and not necessarily invertible mass matrix, +# $u$ the vector of unknowns and $f$ the right-hand-side (RHS). For us $f$ can be interpreted as +# the spatial discretization of all linear and nonlinear operators depending on $u$ and $t$, +# but not on the time derivative of $u$. +# +# ## Some Theory on the Incompressible Navier-Stokes Equations +# +# ### Problem Description in Strong Form +# +# The incompressible Navier-Stokes equations can be stated as the system +# ```math +# \begin{aligned} +# \partial_t v &= \underbrace{\nu \Delta v}_{\text{viscosity}} - \underbrace{(v \cdot \nabla) v}_{\text{advection}} - \underbrace{\nabla p}_{\text{pressure}} \\ +# 0 &= \underbrace{\nabla \cdot v}_{\text{incompressibility}} +# \end{aligned} +# ``` +# where $v$ is the unknown velocity field, $p$ the unknown pressure field, +# $\nu$ the dynamic viscosity and $\Delta$ the Laplacian. In the derivation we assumed +# a constant density of 1 for the fluid and negligible coupling between the velocity components. +# Finally we see that the pressure term appears only in combination with the gradient +# operator, so for any solution $p$ the function $p + c$ is also an admissible solution, if +# we do not impose Dirichlet conditions on the pressure. To resolve this we introduce the +# implicit constraint that $ \int_\Omega p = 0 $. +# +# Our setup is derived from [Turek's DFG benchmark](http://www.mathematik.tu-dortmund.de/~featflow/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark1_re20.html). +# We model a channel with size $0.41 \times 2.2$ and a hole of radius $0.05$ centered at $(0.2, 0.2)$. +# The left side has a parabolic inflow profile, which is ramped up over time, modeled as the time dependent +# Dirichlet condition +# ```math +# v(x,y,t) +# = +# \begin{bmatrix} +# 4 v_{in}(t) y (0.41-y)/0.41^2 \\ +# 0 +# \end{bmatrix} +# ``` +# where $v_{in}(t) = \text{clamp}(t, 0.0, 1.0)$. With a dynamic viscosity of $\nu = 0.001$ +# this is enough to induce turbulence behind the cylinder which leads to vortex shedding. The top and bottom of our +# channel have no-slip conditions, i.e. $v = [0,0]^{\textrm{T}}$, while the right boundary has the do-nothing boundary condtion +# $\nu \partial_{\textrm{n}} v - p n = 0$ to model outflow. With these boundary conditions we can choose the zero solution as a +# feasible initial condition. +# +# ### Derivation of Semi-Discrete Weak Form +# +# By multiplying test functions $\varphi$ and $\psi$ from a suitable test function space on the strong form, +# followed by integrating over the domain and applying partial integration to the pressure and viscosity terms +# we can obtain the following weak form +# ```math +# \begin{aligned} +# \int_\Omega \partial_t v \cdot \varphi &= - \int_\Omega \nu \nabla v : \nabla \varphi - \int_\Omega (v \cdot \nabla) v \cdot \varphi + \int_\Omega p (\nabla \cdot \varphi) + \int_{\partial \Omega_{N}} \underbrace{(\nu \partial_n v - p n )}_{=0} \cdot \varphi \\ +# 0 &= \int_\Omega (\nabla \cdot v) \psi +# \end{aligned} +# ``` +# for all possible test functions from the suitable space. +# +# Now we can discretize the problem as usual with the finite element method +# utilizing Taylor-Hood elements (Q2Q1) to yield a stable discretization in +# mass matrix form: +# ```math +# \underbrace{\begin{bmatrix} +# M_v & 0 \\ +# 0 & 0 +# \end{bmatrix}}_{:=M} +# \begin{bmatrix} +# \mathrm{d}_t\hat{v} \\ +# \mathrm{d}_t\hat{p} +# \end{bmatrix} +# = +# \underbrace{\begin{bmatrix} +# A & B^{\textrm{T}} \\ +# B & 0 +# \end{bmatrix}}_{:=K} +# \begin{bmatrix} +# \hat{v} \\ +# \hat{p} +# \end{bmatrix} +# + +# \begin{bmatrix} +# N(\hat{v}, \hat{v}, \hat{\varphi}) \\ +# 0 +# \end{bmatrix} +# ``` +# Here $M$ is the singular block mass matrix, $K$ is the discretized Stokes operator and $N$ the nonlinear advection term, which +# is also called trilinear form. $\hat{v}$ and $\hat{p}$ represent the time-dependent vectors of nodal values of the discretizations +# of $v$ and $p$ respectively, while $\hat{\varphi}$ is the choice for the test function in the discretization. The hats are dropped +# in the implementation and only stated for clarity in this section. +# +# +# ## Commented Implementation +# +# Now we solve the problem with Ferrite and [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl). What follows is a program spliced with comments. +# The full program, without comments, can be found in the next [section](@ref ns_vs_diffeq-plain-program). +# +# First we load Ferrite and some other packages we need +using Ferrite, SparseArrays, BlockArrays, LinearAlgebra, UnPack +# Since we do not need the complete DifferentialEquations suite, we just load the required ODE infrastructure, which can also handle +# DAEs in mass matrix form. +using OrdinaryDiffEq + +# We start off by defining our only material parameter. +ν = 1.0/1000.0; #dynamic viscosity + +# Next a fine 2D rectangular grid has to be generated. We leave the cell size parametric for flexibility when +# playing around with the code. Note that the mesh is pretty fine, leading to a high memory consumption when +# feeding the equation system to direct solvers. +dim = 2 +cell_scale_factor = 2.0 +x_cells = round(Int, cell_scale_factor*220) +y_cells = round(Int, cell_scale_factor*41) +# CI chokes if the grid is too fine. :) #src +x_cells = round(Int, 55/3) #hide +y_cells = round(Int, 41/3) #hide +grid = generate_grid(Quadrilateral, (x_cells, y_cells), Vec{2}((0.0, 0.0)), Vec{2}((2.2, 0.41))); + +# Next we carve a hole $B_{0.05}(0.2,0.2)$ in the mesh by deleting the cells and update the boundary face sets. +# This code will be replaced once a proper mesh interface is avaliable. +cell_indices = filter(ci->norm(mean(map(i->grid.nodes[i].x-[0.2,0.2], Ferrite.vertices(grid.cells[ci]))))>0.05, 1:length(grid.cells)) +hole_cell_indices = filter(ci->norm(mean(map(i->grid.nodes[i].x-[0.2,0.2], Ferrite.vertices(grid.cells[ci]))))<=0.05, 1:length(grid.cells)); +hole_face_ring = Set{FaceIndex}() +for hci ∈ hole_cell_indices + push!(hole_face_ring, FaceIndex((hci+1, 4))) + push!(hole_face_ring, FaceIndex((hci-1, 2))) + push!(hole_face_ring, FaceIndex((hci-x_cells, 3))) + push!(hole_face_ring, FaceIndex((hci+x_cells, 1))) +end +grid.facesets["hole"] = Set(filter(x->x.idx[1] ∉ hole_cell_indices, collect(hole_face_ring))); +cell_indices_map = map(ci->norm(mean(map(i->grid.nodes[i].x-[0.2,0.2], Ferrite.vertices(grid.cells[ci]))))>0.05 ? indexin([ci], cell_indices)[1] : 0, 1:length(grid.cells)) +grid.cells = grid.cells[cell_indices] +for facesetname in keys(grid.facesets) + grid.facesets[facesetname] = Set(map(fi -> FaceIndex( cell_indices_map[fi.idx[1]] ,fi.idx[2]), collect(grid.facesets[facesetname]))) +end; + +# We test against full development of the flow - so regenerate the grid #src +grid = generate_grid(Quadrilateral, (x_cells, y_cells), Vec{2}((0.0, 0.0)), Vec{2}((0.55, 0.41))); #hide + +# ### Function Space +# To ensure stability we utilize the Taylor-Hood element pair Q2-Q1. +# We have to utilize the same quadrature rule for the pressure as for the velocity, because in the weak form the +# linear pressure term is tested against a quadratic function. +ip_v = Lagrange{dim, RefCube, 2}() +ip_geom = Lagrange{dim, RefCube, 1}() +qr = QuadratureRule{dim, RefCube}(4) +cellvalues_v = CellVectorValues(qr, ip_v, ip_geom); + +ip_p = Lagrange{dim, RefCube, 1}() +cellvalues_p = CellScalarValues(qr, ip_p, ip_geom); + +dh = DofHandler(grid) +push!(dh, :v, dim, ip_v) +push!(dh, :p, 1, ip_p) +close!(dh); + +# ### Boundary Conditions +# As in the DFG benchmark we apply no-slip conditions to the top, bottom and +# cylinder boundary. The no-slip condition states that the velocity of the +# fluid on this portion of the boundary is fixed to be zero. +ch = ConstraintHandler(dh); + +nosplip_face_names = ["top", "bottom", "hole"]; +# No hole for the test present #src +nosplip_face_names = ["top", "bottom"] #hide +∂Ω_noslip = union(getfaceset.((grid, ), nosplip_face_names)...); +noslip_bc = Dirichlet(:v, ∂Ω_noslip, (x, t) -> [0,0], [1,2]) +add!(ch, noslip_bc); + +# The left boundary has a parabolic inflow with peak velocity of 1.0. This +# ensures that for the given geometry the Reynolds number is 100, which +# is already enough to obtain some simple vortex streets. By increasing the +# velocity further we can obtain stronger vortices - which may need additional +# refinement of the grid. +∂Ω_inflow = getfaceset(grid, "left"); + +vᵢₙ(t) = clamp(t, 0.0, 1.0)*1.0 #inflow velocity +vᵢₙ(t) = clamp(t, 0.0, 1.0)*0.3 #hide +parabolic_inflow_profile((x,y),t) = [4*vᵢₙ(t)*y*(0.41-y)/0.41^2,0] +inflow_bc = Dirichlet(:v, ∂Ω_inflow, parabolic_inflow_profile, [1,2]) +add!(ch, inflow_bc); + +# The outflow boundary condition has been applied on the right side of the +# cylinder when the weak form has been derived by setting the boundary integral +# to zero. It is also called the do-nothing condition. Other outflow conditions +# are also possible. +∂Ω_free = getfaceset(grid, "right"); + +close!(ch) +update!(ch, 0.0); + +# ### Linear System Assembly +# Next we describe how the block mass matrix and the Stokes matrix are assembled. +# +# For the block mass matrix $M$ we remember that only the first equation had a time derivative +# and that the block mass matrix corresponds to the term arising from discretizing the time +# derivatives. Hence, only the upper left block has non-zero components. +function assemble_mass_matrix(cellvalues_v::CellVectorValues{dim}, cellvalues_p::CellScalarValues{dim}, M::SparseMatrixCSC, dh::DofHandler) where {dim} + ## Allocate a buffer for the local matrix and some helpers, together with the assembler. + n_basefuncs_v = getnbasefunctions(cellvalues_v) + n_basefuncs_p = getnbasefunctions(cellvalues_p) + n_basefuncs = n_basefuncs_v + n_basefuncs_p + v▄, p▄ = 1, 2 + Mₑ = PseudoBlockArray(zeros(n_basefuncs, n_basefuncs), [n_basefuncs_v, n_basefuncs_p], [n_basefuncs_v, n_basefuncs_p]) + + ## It follows the assembly loop as explained in the basic tutorials. + mass_assembler = start_assemble(M) + @inbounds for cell in CellIterator(dh) + fill!(Mₑ, 0) + Ferrite.reinit!(cellvalues_v, cell) + + for q_point in 1:getnquadpoints(cellvalues_v) + dΩ = getdetJdV(cellvalues_v, q_point) + ## Remember that we assemble a vector mass term, hence the dot product. + for i in 1:n_basefuncs_v + φᵢ = shape_value(cellvalues_v, q_point, i) + for j in 1:n_basefuncs_v + φⱼ = shape_value(cellvalues_v, q_point, j) + Mₑ[BlockIndex((v▄, v▄), (i, j))] += φᵢ ⋅ φⱼ * dΩ + end + end + end + assemble!(mass_assembler, celldofs(cell), Mₑ) + end + + return M +end; + +# Next we discuss the assembly of the Stokes matrix. +# Remember that we use the same function spaces for trial and test, hence the +# matrix has the following block form +# ```math +# K = \begin{bmatrix} +# A & B^{\textrm{T}} \\ +# B & 0 +# \end{bmatrix} +# ``` +# which is also called saddle point matrix. These problems are known to have +# a non-trivial kernel, which is a reflection of the strong form as discussed +# in the theory portion if this example. +function assemble_stokes_matrix(cellvalues_v::CellVectorValues{dim}, cellvalues_p::CellScalarValues{dim}, ν, K::SparseMatrixCSC, dh::DofHandler) where {dim} + ## Again, some buffers and helpers + n_basefuncs_v = getnbasefunctions(cellvalues_v) + n_basefuncs_p = getnbasefunctions(cellvalues_p) + n_basefuncs = n_basefuncs_v + n_basefuncs_p + v▄, p▄ = 1, 2 + Kₑ = PseudoBlockArray(zeros(n_basefuncs, n_basefuncs), [n_basefuncs_v, n_basefuncs_p], [n_basefuncs_v, n_basefuncs_p]) + + ## Assembly loop + stiffness_assembler = start_assemble(K) + @inbounds for cell in CellIterator(dh) + ## Don't forget to initialize everything + fill!(Kₑ, 0) + + Ferrite.reinit!(cellvalues_v, cell) + Ferrite.reinit!(cellvalues_p, cell) + + for q_point in 1:getnquadpoints(cellvalues_v) + dΩ = getdetJdV(cellvalues_v, q_point) + # Assemble local viscosity block of $A$ + #+ + for i in 1:n_basefuncs_v + ∇φᵢ = shape_gradient(cellvalues_v, q_point, i) + for j in 1:n_basefuncs_v + ∇φⱼ = shape_gradient(cellvalues_v, q_point, j) + Kₑ[BlockIndex((v▄, v▄), (i, j))] -= ν * ∇φᵢ ⊡ ∇φⱼ * dΩ + end + end + # Assemble local pressure and incompressibility blocks of $B^{\textrm{T}}$ and $B$. + #+ + for j in 1:n_basefuncs_p + ψ = shape_value(cellvalues_p, q_point, j) + for i in 1:n_basefuncs_v + divφ = shape_divergence(cellvalues_v, q_point, i) + Kₑ[BlockIndex((v▄, p▄), (i, j))] += (divφ * ψ) * dΩ + Kₑ[BlockIndex((p▄, v▄), (j, i))] += (ψ * divφ) * dΩ + end + end + end + + ## Assemble `Kₑ` into the Stokes matrix `K`. + assemble!(stiffness_assembler, celldofs(cell), Kₑ) + end + return K +end; + +# ### Solution of the semi-discretized system via DifferentialEquations.jl +# First we assemble the linear portions for efficiency. These matrices are +# assumed to be constant over time. +T = 10.0 +Δt₀ = 0.01 +Δt_save = 0.1 + +M = create_sparsity_pattern(dh); +M = assemble_mass_matrix(cellvalues_v, cellvalues_p, M, dh); + +K = create_sparsity_pattern(dh); +K = assemble_stokes_matrix(cellvalues_v, cellvalues_p, ν, K, dh); + +# These are our initial conditions. We start from the zero solution, because it +# is trivially admissible if the Dirichlet conditions are zero everywhere on the +# Dirichlet boundary for $t=0$. Note that the time stepper is also doing fine if the +# Dirichlet condition is non-zero and not too pathological. +u₀ = zeros(ndofs(dh)) +apply!(u₀, ch); + +# DifferentialEquations assumes dense matrices by default, which is not +# feasible for semi-discretization of finite element models. We communicate +# that a sparse matrix with specified pattern should be utilized through the +# `jac_prototyp` argument. It is simple to see that the Jacobian and the +# stiffness matrix share the same sparsity pattern, since they share the +# same relation between trial and test functions. +jac_sparsity = sparse(K); + +# To apply the nonlinear portion of the Navier-Stokes problem we simply hand +# over the dof handler and cell values to the right-hand-side (RHS) as a parameter. +# Further the pre-assembled linear part (which is time independent) is +# passed to save some runtime. To apply the time-dependent Dirichlet BCs, we +# also hand over the constraint handler. +# The basic idea to apply the Dirichlet BCs consistently is that we copy the +# current solution `u`, apply the Dirichlet BCs on the copy, evaluate the +# discretized RHS of the Navier-Stokes equations with this vector +# and finally set the RHS to zero on every constraint. This way we obtain a +# correct solution for all dofs which are not Dirichlet constrained. These +# dofs are then corrected in a post-processing step, when evaluating the +# solution vector at specific time points. +# It should be finally noted that this **trick does not work** out of the box +# **for constraining algebraic portion** of the DAE, i.e. if we would like to +# put a Dirichlet BC on pressure dofs. As a workaround we have to set $f_{\textrm{i}} = 1$ +# instead of $f_{\textrm{i}} = 0$, because otherwise the equation system gets singular. +# This is obvious when we remember that our mass matrix is zero for these +# dofs, such that we obtain the equation $0 \cdot \mathrm{d}_t p_{\textrm{i}} = 1 \cdot p_{\textrm{i}}$, which +# now has a unique solution. +struct RHSparams + K::SparseMatrixCSC + ch::ConstraintHandler + dh::DofHandler + cellvalues_v::CellVectorValues +end +p = RHSparams(K, ch, dh, cellvalues_v) + +function navierstokes!(du,u_uc,p,t) + # Unpack the struct to save some allocations. + #+ + @unpack K,ch,dh,cellvalues_v = p + + # We start by applying the time-dependent Dirichlet BCs. Note that we are + # not allowed to mutate `u_uc`! We also can not pre-allocate this variable + # if we want to use AD to derive the Jacobian matrix, which appears in the + # utilized implicit Euler. If we hand over the Jacobian analytically to + # the solver, or when utilizing a method which does not require building the + # Jacobian, then we could also hand over a buffer for `u` in our RHSparams + # structure to save the allocations made here. + #+ + u = u_uc + update!(ch, t) + apply!(u, ch) + + # Now we apply the rhs of the Navier-Stokes equations + #+ + ## Linear contribution (Stokes operator) + du .= K * u + + ## nonlinear contribution + n_basefuncs = getnbasefunctions(cellvalues_v) + for cell in CellIterator(dh) + Ferrite.reinit!(cellvalues_v, cell) + all_celldofs = celldofs(cell) + v_celldofs = all_celldofs[dof_range(dh, :v)] + v_cell = u[v_celldofs] + for q_point in 1:getnquadpoints(cellvalues_v) + dΩ = getdetJdV(cellvalues_v, q_point) + ∇v = function_gradient(cellvalues_v, q_point, v_cell) + v = function_value(cellvalues_v, q_point, v_cell) + for j in 1:n_basefuncs + φⱼ = shape_value(cellvalues_v, q_point, j) + # Note that in Tensors.jl the definition $\textrm{grad} v = \nabla v$ holds. + # With this information it can be quickly shown in index notation that + # ```math + # [(v \cdot \nabla) v]_{\textrm{i}} = v_{\textrm{j}} (\partial_{\textrm{j}} v_{\textrm{i}}) = [v (\nabla v)^{\textrm{T}}]_{\textrm{i}} + # ``` + # where we should pay attentation to the transpose of the gradient. + #+ + du[v_celldofs[j]] -= v ⋅ ∇v' ⋅ φⱼ * dΩ + end + end + end + + # For now we have to ingore the evolution of the Dirichlet BCs. + # The DBC dofs in the solution vector will be corrected in a post-processing step. + #+ + apply_zero!(du, ch) +end; +# Finally, together with our pre-assembled mass matrix, we are now able to +# define our problem in mass matrix form. +rhs = ODEFunction(navierstokes!, mass_matrix=M; jac_prototype=jac_sparsity) +problem = ODEProblem(rhs, u₀, (0.0,T), p); + +# Now we can put everything together by specifying how to solve the problem. +# We want to use the adaptive implicit Euler method with our custom linear +# solver, which helps in the enforcement of the Dirichlet BCs. Further we +# enable the progress bar with the `progess` and `progress_steps` arguments. +# Finally we have to communicate the time step length and initialization +# algorithm. Since we start with a valid initial state we do not use one of +# DifferentialEquations.jl initialization algorithms. +# NOTE: At the time of writing this [no Hessenberg index 2 initialization is implemented](https://github.com/SciML/OrdinaryDiffEq.jl/issues/1019). +# +# To visualize the result we export the grid and our fields +# to VTK-files, which can be viewed in [ParaView](https://www.paraview.org/) +# by utilizing the corresponding pvd file. +timestepper = ImplicitEuler() +integrator = init( + problem, timestepper, initializealg=NoInit(), dt=Δt₀, + adaptive=true, abstol=1e-3, reltol=1e-3, + progress=true, progress_steps=1, + saveat=Δt_save); + +pvd = paraview_collection("vortex-street.pvd"); +integrator = TimeChoiceIterator(integrator, 0.0:Δt_save:T) +for (u_uc,t) in integrator + # We ignored the Dirichlet constraints in the solution vector up to now, + # so we have to bring them back now. + #+ + update!(ch, t) + u = u_uc + apply!(u, ch) + #compress=false flag because otherwise each vtk file will be stored in memory + vtk_grid("vortex-street-$t.vtu", dh; compress=false) do vtk + vtk_point_data(vtk,dh,u) + vtk_save(vtk) + pvd[t] = vtk + end +end +vtk_save(pvd); + +# Test the result for full proper development of the flow #src +using Test #hide +function compute_divergence(dh, u, cellvalues_v) #hide + divv = 0.0 #hide + @inbounds for (i,cell) in enumerate(CellIterator(dh)) #hide + Ferrite.reinit!(cellvalues_v, cell) #hide + for q_point in 1:getnquadpoints(cellvalues_v) #hide + dΩ = getdetJdV(cellvalues_v, q_point) #hide + #hide + all_celldofs = celldofs(cell) #hide + v_celldofs = all_celldofs[dof_range(dh, :v)] #hide + v_cell = u[v_celldofs] #hide + #hide + divv += function_divergence(cellvalues_v, q_point, v_cell) * dΩ #hide + end #hide + end #hide + return divv #hide +end #hide +@testset "INS OrdinaryDiffEq" begin #hide + u = integrator.integrator.u #hide + Δdivv = abs(compute_divergence(dh, u, cellvalues_v)) #hide + @test isapprox(Δdivv, 0.0, atol=1e-12) #hide + #hide + Δv = 0.0 #hide + for cell in CellIterator(dh) #hide + Ferrite.reinit!(cellvalues_v, cell) #hide + all_celldofs = celldofs(cell) #hide + v_celldofs = all_celldofs[dof_range(dh, :v)] #hide + v_cell = u[v_celldofs] #hide + coords = getcoordinates(cell) #hide + for q_point in 1:getnquadpoints(cellvalues_v) #hide + dΩ = getdetJdV(cellvalues_v, q_point) #hide + coords_qp = spatial_coordinate(cellvalues_v, q_point, coords) #hide + v = function_value(cellvalues_v, q_point, v_cell) #hide + Δv += norm(v - parabolic_inflow_profile(coords_qp, T))^2*dΩ #hide + end #hide + end #hide + @test isapprox(sqrt(Δv), 0.0, atol=1e-3) #hide +end; #hide + +#md # ## [Plain Program](@id ns_vs_diffeq-plain-program) +#md # +#md # Below follows a version of the program without any comments. +#md # The file is also available here: [ns_vs_diffeq.jl](ns_vs_diffeq.jl) +#md # +#md # ```julia +#md # @__CODE__ +#md # ``` \ No newline at end of file diff --git a/docs/src/literate/plasticity.jl b/docs/src/literate/plasticity.jl index ddd462d622..70b1edac24 100644 --- a/docs/src/literate/plasticity.jl +++ b/docs/src/literate/plasticity.jl @@ -33,7 +33,7 @@ # ### Material parameters and state variables # # Start by loading some necessary packages -using Ferrite, SparseArrays, LinearAlgebra, Printf +using Ferrite, Tensors, SparseArrays, LinearAlgebra, Printf # We define a J₂-plasticity-material, containing material parameters and the elastic # stiffness Dᵉ (since it is constant) @@ -62,37 +62,22 @@ end; #md # material parameters ``E`` and ``ν`` - simply as a convenience for the user. #md # -# Define a `struct` to store the material state. -mutable struct MaterialState{T, S <: SecondOrderTensor{3, T}} +# Define a `struct` to store the material state for a Gauss point. +struct MaterialState{T, S <: SecondOrderTensor{3, T}} ## Store "converged" values ϵᵖ::S # plastic strain σ::S # stress k::T # hardening variable - - ## Store temporary values used during equilibrium iterations - temp_ϵᵖ::S - temp_σ::S - temp_k::T end # Constructor for initializing a material state. Every quantity is set to zero. function MaterialState() return MaterialState( - zero(SymmetricTensor{2, 3}), - zero(SymmetricTensor{2, 3}), - 0.0, zero(SymmetricTensor{2, 3}), zero(SymmetricTensor{2, 3}), 0.0) end -# Next, we define a method to update the material state after equilibrium has -# been found. This will be called at the end of each time-step. -function update_state!(state::MaterialState) - state.ϵᵖ = state.temp_ϵᵖ - state.σ = state.temp_σ - state.k = state.temp_k -end; # For later use, during the post-processing step, we define a function to # compute the von Mises effective stress. function vonMises(σ) @@ -104,25 +89,23 @@ end; # # This is the actual method which computes the stress and material tangent # stiffness in a given integration point. -# Input is the current strain and material state. +# Input is the current strain and the material state from the previous timestep. function compute_stress_tangent(ϵ::SymmetricTensor{2, 3}, material::J2Plasticity, state::MaterialState) ## unpack some material parameters G = material.G - K = material.K H = material.H ## We use (•)ᵗ to denote *trial*-values σᵗ = material.Dᵉ ⊡ (ϵ - state.ϵᵖ) # trial-stress sᵗ = dev(σᵗ) # deviatoric part of trial-stress J₂ = 0.5 * sᵗ ⊡ sᵗ # second invariant of sᵗ - σᵗₑ = sqrt(3.0*J₂) # effetive trial-stress (von Mises stress) + σᵗₑ = sqrt(3.0*J₂) # effective trial-stress (von Mises stress) σʸ = material.σ₀ + H * state.k # Previous yield limit φᵗ = σᵗₑ - σʸ # Trial-value of the yield surface if φᵗ < 0.0 # elastic loading - state.temp_σ = σᵗ - return state.temp_σ, material.Dᵉ + return σᵗ, material.Dᵉ, MaterialState(state.ϵᵖ, σᵗ, state.k) else # plastic loading h = H + 3G μ = φᵗ / h # plastic multiplier @@ -143,12 +126,11 @@ function compute_stress_tangent(ϵ::SymmetricTensor{2, 3}, material::J2Plasticit Dtemp(i,j,k,l) = -2G*b * Q(i,j,k,l) - 9G^2 / (h*σₑ^2) * s[i,j]*s[k,l] D = material.Dᵉ + SymmetricTensor{4, 3}(Dtemp) - ## Store outputs in the material state - Δϵᵖ = 3/2 *μ / σₑ*s # plastic strain - state.temp_ϵᵖ = state.ϵᵖ + Δϵᵖ # plastic strain - state.temp_k = state.k + μ # hardening variable - state.temp_σ = σ # updated stress - return state.temp_σ, D + ## Return new state + Δϵᵖ = 3/2 * μ / σₑ * s # plastic strain + ϵᵖ = state.ϵᵖ + Δϵᵖ # plastic strain + k = state.k + μ # hardening variable + return σ, D, MaterialState(ϵᵖ, σ, k) end end @@ -196,20 +178,22 @@ end; # * Tangent stiffness `K` function doassemble(cellvalues::CellVectorValues{dim}, facevalues::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid, - dh::DofHandler, material::J2Plasticity, u, states, t) where {dim} + dh::DofHandler, material::J2Plasticity, u, states, states_old, t) where {dim} r = zeros(ndofs(dh)) assembler = start_assemble(K, r) nu = getnbasefunctions(cellvalues) re = zeros(nu) # element residual vector ke = zeros(nu, nu) # element tangent matrix - for (cell, state) in zip(CellIterator(dh), states) + for (i, cell) in enumerate(CellIterator(dh)) fill!(ke, 0) fill!(re, 0) eldofs = celldofs(cell) ue = u[eldofs] + state = @view states[:, i] + state_old = @view states_old[:, i] assemble_cell!(ke, re, cell, cellvalues, facevalues, grid, material, - ue, state, t) + ue, state, state_old, t) assemble!(assembler, eldofs, re, ke) end return K, r @@ -221,23 +205,21 @@ end #md # and then symmetrize it. #md # function assemble_cell!(Ke, re, cell, cellvalues, facevalues, grid, material, - ue, state, t) + ue, state, state_old, t) n_basefuncs = getnbasefunctions(cellvalues) reinit!(cellvalues, cell) for q_point in 1:getnquadpoints(cellvalues) ## For each integration point, compute stress and material stiffness - ∇u = function_gradient(cellvalues, q_point, ue) - ϵ = symmetric(∇u) # Total strain - σ, D = compute_stress_tangent(ϵ, material, state[q_point]) + ϵ = function_symmetric_gradient(cellvalues, q_point, ue) # Total strain + σ, D, state[q_point] = compute_stress_tangent(ϵ, material, state_old[q_point]) dΩ = getdetJdV(cellvalues, q_point) for i in 1:n_basefuncs - δϵ = symmetric(shape_gradient(cellvalues, q_point, i)) - + δϵ = shape_symmetric_gradient(cellvalues, q_point, i) re[i] += (δϵ ⊡ σ) * dΩ # add internal force to residual - for j in 1:i - Δϵ = symmetric(shape_gradient(cellvalues, q_point, j)) + for j in 1:i # loop only over lower half + Δϵ = shape_symmetric_gradient(cellvalues, q_point, j) Ke[i, j] += δϵ ⊡ D ⊡ Δϵ * dΩ end end @@ -308,10 +290,8 @@ function solve() ## Create material states. One array for each cell, where each element is an array of material- ## states - one for each integration point nqp = getnquadpoints(cellvalues) - states = [[MaterialState() for _ in 1:nqp] for _ in 1:getncells(grid)] - - ## states = [MaterialState() for _ in 1:nqp for _ in 1:getncells(grid)] - ## temp_states = [MaterialState() for _ in 1:nqp for _ in 1:getncells(grid)] + states = [MaterialState() for _ in 1:nqp, _ in 1:getncells(grid)] + states_old = [MaterialState() for _ in 1:nqp, _ in 1:getncells(grid)] ## Newton-Raphson loop NEWTON_TOL = 1 # 1 N @@ -332,7 +312,7 @@ function solve() break end K, r = doassemble(cellvalues, facevalues, K, grid, dh, material, u, - states, traction); + states, states_old, traction); norm_r = norm(r[Ferrite.free_dofs(dbcs)]) print("Iteration: $newton_itr \tresidual: $(@sprintf("%.8f", norm_r))\n") @@ -345,10 +325,9 @@ function solve() u -= Δu end - ## Update all the material states after we have reached equilibrium - for cell_states in states - foreach(update_state!, cell_states) - end + ## Update the old states with the converged values for next timestep + states_old .= states + u_max[timestep] = max(abs.(u)...) # maximum displacement in current timestep end @@ -358,7 +337,7 @@ function solve() ## The following is a quick (and dirty) way of extracting average cell data for export. mises_values = zeros(getncells(grid)) κ_values = zeros(getncells(grid)) - for (el, cell_states) in enumerate(states) + for (el, cell_states) in enumerate(eachcol(states)) for state in cell_states mises_values[el] += vonMises(state.σ) κ_values[el] += state.k*material.H @@ -386,7 +365,7 @@ plot( vcat(0.0, traction_magnitude), linewidth=2, title="Traction-displacement", - label=[""], + label=nothing, markershape=:auto ) ylabel!("Traction [Pa]") @@ -398,7 +377,7 @@ xlabel!("Maximum deflection [m]") ## test the result #src using Test #src -@test norm(u_max[end]) ≈ 0.2544526451 #src +@test norm(u_max[end]) ≈ 0.254452645 #src #md # ## [Raw source](@id plasticity-raw-code) #md # diff --git a/docs/src/literate/postprocessing.jl b/docs/src/literate/postprocessing.jl new file mode 100644 index 0000000000..76c51124cf --- /dev/null +++ b/docs/src/literate/postprocessing.jl @@ -0,0 +1,126 @@ +# # Postprocessing +# +# ![](heat_square_fluxes.png) + + +# ## Introduction +# +# After running a simulation, we usually want to visualize the results in different ways. +# The `L2Projector` and the `PointEvalHandler` build a pipeline for doing so. With the `L2Projector`, +# integration point quantities can be projected to the nodes. The `PointEvalHandler` enables evaluation of +# the finite element approximated function in any coordinate in the domain. Thus with the combination of both functionalities, +# both nodal quantities and integration point quantities can be evaluated in any coordinate, allowing for example +# cut-planes through 3D structures or cut-lines through 2D-structures. +# +# This example continues from the Heat equation example, where the temperature field was +# determined on a square domain. In this example, we first compute the heat flux in each +# integration point (based on the solved temperature field) and then we do an L2-projection +# of the fluxes to the nodes of the mesh. By doing this, we can more easily visualize +# integration points quantities. Finally, we visualize the temperature field and the heat fluxes along a cut-line. +# +# The L2-projection is defined as follows: Find projection ``q(\boldsymbol{x}) \in L_2(\Omega)`` such that +# ```math +# \int v q \ \mathrm{d}\Omega = \int v d \ \mathrm{d}\Omega \quad \forall v \in L_2(\Omega), +# ``` +# where ``d`` is the quadrature data to project. Since the flux is a vector the projection function +# will be solved with multiple right hand sides, e.g. with ``d = q_x`` and ``d = q_y`` for this 2D problem. +# +# Ferrite has functionality for doing much of this automatically, as displayed in the code below. +# In particular [`L2Projector`](@ref) for assembling the left hand side, and +# [`project`](@ref) for assembling the right hand sides and solving for the projection. + +# ## Implementation +# +# Start by simply running the Heat equation example to solve the problem +include("heat_equation.jl"); + + +# Next we define a function that computes the heat flux for each integration point in the domain. +# Fourier's law is adopted, where the conductivity tensor is assumed to be isotropic with unit +# conductivity ``\lambda = 1 ⇒ q = - \nabla u``, where ``u`` is the temperature. +function compute_heat_fluxes(cellvalues::CellScalarValues{dim,T}, dh::DofHandler, a) where {dim,T} + + n = getnbasefunctions(cellvalues) + cell_dofs = zeros(Int, n) + nqp = getnquadpoints(cellvalues) + + ## Allocate storage for the fluxes to store + q = [Vec{2,T}[] for _ in 1:getncells(dh.grid)] + + for (cell_num, cell) in enumerate(CellIterator(dh)) + q_cell = q[cell_num] + celldofs!(cell_dofs, dh, cell_num) + aᵉ = a[cell_dofs] + reinit!(cellvalues, cell) + + for q_point in 1:nqp + q_qp = - function_gradient(cellvalues, q_point, aᵉ) + push!(q_cell, q_qp) + end + end + return q +end +#md nothing # hide + +# Now call the function to get all the fluxes. +q_gp = compute_heat_fluxes(cellvalues, dh, u); + +# Next, create an [`L2Projector`](@ref) using the same interpolation as was used to approximate the +# temperature field. On instantiation, the projector assembles the coefficient matrix `M` and +# computes the Cholesky factorization of it. By doing so, the projector can be reused without +# having to invert `M` every time. +projector = L2Projector(ip, grid); + +# Project the integration point values to the nodal values +q_projected = project(projector, q_gp, qr; project_to_nodes=false); # TODO: this should be default. + + +# ## Exporting to VTK +# To visualize the heat flux, we export the projected field `q_projected` +# to a VTK-file, which can be viewed in e.g. [ParaView](https://www.paraview.org/). +vtk_grid("heat_equation_flux", grid) do vtk + ## TODO: This doesn't work (correctly) yet (https://github.com/Ferrite-FEM/Ferrite.jl/issues/278) + vtk_point_data(vtk, q_projected, "q") +end; + +# ## Point Evaluation +# ![](heat_square_pointevaluation.png) + +# Consider a cut-line through the domain, like the black line in the figure above. +# We will evaluate the temperature and the heat flux distribution along a horizontal line. +points = [Vec((x, 0.75)) for x in range(-1.0, 1.0, length=101)]; + +# First, we need to generate a `PointEvalHandler`. This will find and store the cells +# containing the input points. +ph = PointEvalHandler(grid, points); + +# After the L2-Projection, the heat fluxes `q_projected` are stored in the DoF-ordering +# determined by the projector's internal DoFHandler, so to evalute the flux `q` at our points +# we give the `PointEvalHandler`, the `L2Projector` and the values `q_projected`. +q_points = get_point_values(ph, projector, q_projected); + +# We can also extract the field values, here the temperature, right away from the result +# vector of the simulation, that is stored in `u`. These values are stored in the order of +# our initial DofHandler so the input is not the `PointEvalHandler`, the original `DofHandler`, +# the dof-vector `u`, and (optionally for single-field problems) the name of the field. +# from the `L2Projection`, the values are stored in the order of the degrees of freedom. +u_points = Ferrite.get_point_values(ph, dh, u, :u); + +# Now, we can plot the temperature and flux values with the help of any plotting library, e.g. Plots.jl. +# To do this, we need to import the package: +import Plots + +# Firstly, we are going to plot the temperature values along the given line. +Plots.plot(getindex.(points,1), u_points, xlabel="x (coordinate)", ylabel="u (temperature)", label=nothing) + +# Secondly, the horizontal heat flux (i.e. the first component of the heat flux vector) is plotted. +Plots.plot(getindex.(points,1), getindex.(q_points,1), xlabel="x (coordinate)", ylabel="q_x (flux in x-direction)", label=nothing) + +#md # ## [Plain Program](@id postprocessing-plain-program) +#md # +#md # Below follows a version of the program without any comments. +#md # The file is also available here: [postprocessing.jl](postprocessing.jl) +#md # +#md # ```julia +#md # @__CODE__ +#md # ``` diff --git a/docs/src/reference/export.md b/docs/src/reference/export.md index b0557db5c9..908ad9bc89 100644 --- a/docs/src/reference/export.md +++ b/docs/src/reference/export.md @@ -1,14 +1,30 @@ ```@meta DocTestSetup = :(using Ferrite) ``` +# Postprocessing + +## Project to nodes +```@docs +L2Projector +project +``` -# Export +# Postprocessing ```@docs -vtk_grid +PointEvalHandler +get_point_values ``` ```@docs -L2Projector -project +reshape_to_nodes ``` + +## VTK Export + +```@docs +vtk_grid(filename::AbstractString, grid::Grid{dim,C,T}; compress::Bool) where {dim,C,T} +vtk_point_data(vtk::WriteVTK.DatasetFile, data::Union{Vector{SymmetricTensor{2,dim,T,M}}},name::AbstractString) where {dim,T,M} +vtk_point_data(vtk::WriteVTK.DatasetFile, data::Union{ Vector{Tensor{order,dim,T,M}}, Vector{SymmetricTensor{order,dim,T,M}}}, name::AbstractString) where {order,dim,T,M} +vtk_cellset +``` \ No newline at end of file diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index 01d5f2b931..5c948d5612 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -69,6 +69,12 @@ getfieldinterpolation(dh::DofHandler, field_idx::Int) = dh.field_interpolations[ getfielddim(dh::DofHandler, field_idx::Int) = dh.field_dims[field_idx] getbcvalue(dh::DofHandler, field_idx::Int) = dh.bc_values[field_idx] +function getfielddim(dh::DofHandler, name::Symbol) + field_pos = findfirst(i->i == name, getfieldnames(dh)) + field_pos === nothing && error("did not find field $name") + return dh.field_dims[field_pos] +end + """ dof_range(dh:DofHandler, field_name) @@ -387,26 +393,46 @@ function WriteVTK.vtk_grid(filename::AbstractString, dh::AbstractDofHandler; com vtk_grid(filename, dh.grid; compress=compress) end -# Exports the FE field `u` to `vtkfile` -function WriteVTK.vtk_point_data(vtkfile, dh::DofHandler, u::Vector, suffix="") - for f in 1:nfields(dh) - @debug println("exporting field $(dh.field_names[f])") - field_dim = dh.field_dims[f] - space_dim = field_dim == 2 ? 3 : field_dim - data = fill(0.0, space_dim, getnnodes(dh.grid)) - offset = field_offset(dh, dh.field_names[f]) - for cell in CellIterator(dh) - _celldofs = celldofs(cell) - counter = 1 - for node in getnodes(cell) - for d in 1:dh.field_dims[f] - data[d, node] = u[_celldofs[counter + offset]] - @debug println(" exporting $(u[_celldofs[counter + offset]]) for dof#$(_celldofs[counter + offset])") - counter += 1 - end +""" + reshape_to_nodes(dh::AbstractDofHandler, u::Vector{T}, fieldname::Symbol) where T + +Reshape the entries of the dof-vector `u` which correspond to the field `fieldname` in nodal order. +Return a matrix with a column for every node and a row for every dimension of the field. +For superparametric fields only the entries corresponding to nodes of the grid will be returned. Do not use this function for subparametric approximations. +""" +function reshape_to_nodes(dh::DofHandler, u::Vector{T}, fieldname::Symbol) where T + # make sure the field exists + fieldname ∈ Ferrite.getfieldnames(dh) || error("Field $fieldname not found.") + + field_idx = findfirst(i->i==fieldname, getfieldnames(dh)) + offset = field_offset(dh, fieldname) + field_dim = getfielddim(dh, field_idx) + + space_dim = field_dim == 2 ? 3 : field_dim + data = fill(zero(T), space_dim, getnnodes(dh.grid)) + + reshape_field_data!(data, dh, u, offset, field_dim) + + return data +end + +function reshape_field_data!(data::Matrix{T}, dh::AbstractDofHandler, u::Vector{T}, field_offset::Int, field_dim::Int, cellset=Set{Int}(1:getncells(dh.grid))) where T + + _celldofs = Vector{Int}(undef, ndofs_per_cell(dh, first(cellset))) + for cell in CellIterator(dh, collect(cellset)) + celldofs!( _celldofs, cell) + counter = 1 + for node in getnodes(cell) + for d in 1:field_dim + data[d, node] = u[_celldofs[counter + field_offset]] + @debug println(" exporting $(u[_celldofs[counter + field_offset]]) for dof#$(_celldofs[counter + field_offset])") + counter += 1 + end + if field_dim == 2 + # paraview requires 3D-data so pad with zero + data[3, node] = 0 end end - vtk_point_data(vtkfile, data, string(dh.field_names[f], suffix)) end - return vtkfile -end + return data +end \ No newline at end of file diff --git a/src/Dofs/MixedDofHandler.jl b/src/Dofs/MixedDofHandler.jl index f387f34de5..90abc17d19 100644 --- a/src/Dofs/MixedDofHandler.jl +++ b/src/Dofs/MixedDofHandler.jl @@ -406,7 +406,7 @@ function Ferrite.dof_range(fh::FieldHandler, field_name::Symbol) offset = Ferrite.field_offset(fh, field_name) field_interpolation = fh.fields[f].interpolation field_dim = fh.fields[f].dim - n_field_dofs = getnbasefunctions(field_interpolation) * field_dim + n_field_dofs = getnbasefunctions(field_interpolation)::Int * field_dim return (offset+1):(offset+n_field_dofs) end @@ -414,3 +414,23 @@ find_field(dh::MixedDofHandler, field_name::Symbol) = find_field(first(dh.fieldh field_offset(dh::MixedDofHandler, field_name::Symbol) = field_offset(first(dh.fieldhandlers), field_name) getfieldinterpolation(dh::MixedDofHandler, field_idx::Int) = dh.fieldhandlers[1].fields[field_idx].interpolation getfielddim(dh::MixedDofHandler, field_idx::Int) = dh.fieldhandlers[1].fields[field_idx].dim + +function reshape_to_nodes(dh::MixedDofHandler, u::Vector{Float64}, fieldname::Symbol) + + # make sure the field exists + fieldname ∈ Ferrite.getfieldnames(dh) || error("Field $fieldname not found.") + + field_dim = getfielddim(dh, fieldname) + space_dim = field_dim == 2 ? 3 : field_dim + data = fill(NaN, space_dim, getnnodes(dh.grid)) # set default value + + for fh in dh.fieldhandlers + # check if this fh contains this field, otherwise continue to the next + field_pos = findfirst(i->i == fieldname, getfieldnames(fh)) + field_pos === nothing && continue + offset = field_offset(fh, fieldname) + + reshape_field_data!(data, dh, u, offset, field_dim, fh.cellset) + end + return data +end \ No newline at end of file diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl index d049120a15..553ea9be62 100644 --- a/src/Export/VTK.jl +++ b/src/Export/VTK.jl @@ -31,11 +31,43 @@ function WriteVTK.vtk_grid(filename::AbstractString, grid::Grid{dim,C,T}; compre return vtk_grid(filename, coords, cls; compress=compress) end +""" + vtk_point_data(vtk, data::Vector{<:SymmetricTensor{2}}, name) + +Write the second order tensor field `data` to the vtk file. Twodimensional tensors are padded with zeros. +The tensor field is written such that Paraview recognizes the tensor components, where `XX` corresponds to the `[1,1]` component and so on. +""" +function WriteVTK.vtk_point_data( + vtk::WriteVTK.DatasetFile, + data::Union{ + Vector{SymmetricTensor{2,dim,T,M}} + }, + name::AbstractString + ) where {dim,T,M} + + if dim == 1 + noutputs = 1 + sort_paraview = [1] + elseif dim == 2 + noutputs = 6 + sort_paraview = [1,4,2] + elseif dim == 3 + noutputs = 6 + sort_paraview = [1,4,6,2,5,3] + end + + npoints = length(data) + out = zeros(T, noutputs, npoints) + out[sort_paraview, :] .= reshape(reinterpret(T, data), (M, npoints)) + + return vtk_point_data(vtk, out, name) +end + """ vtk_point_data(vtk, data::Vector{<:Tensor}, name) Write the tensor field data to the vtk file. Only writes the tensor values available in `data`. -In the vtu-file, ordering of the tensor components is coulumn-wise (just like Julia). +In the vtu-file, ordering of the tensor components is column-wise (just like Julia). [1 2 3 4] => 1, 3, 2, 4 """ @@ -100,47 +132,14 @@ the cell is in the set and 0 otherwise. vtk_cellset(vtk::WriteVTK.DatasetFile, grid::AbstractGrid, cellset::String) = vtk_cellset(vtk, grid, [cellset]) -import Ferrite.field_offset -function WriteVTK.vtk_point_data(vtkfile, dh::MixedDofHandler, u::Vector, suffix="") +function WriteVTK.vtk_point_data(vtkfile, dh::AbstractDofHandler, u::Vector, suffix="") fieldnames = Ferrite.getfieldnames(dh) # all primary fields for name in fieldnames - @debug println("exporting field $(name)") - field_dim = getfielddim(dh, name) - space_dim = field_dim == 2 ? 3 : field_dim - data = fill(NaN, space_dim, getnnodes(dh.grid)) # set default value - - for fh in dh.fieldhandlers - # check if this fh contains this field, otherwise continue to the next - field_pos = findfirst(i->i == name, getfieldnames(fh)) - field_pos === nothing && continue - - cellnumbers = sort(collect(fh.cellset)) # TODO necessary to have them ordered? - offset = field_offset(fh, name) - - for cellnum in cellnumbers - cell = dh.grid.cells[cellnum] - n = ndofs_per_cell(dh, cellnum) - eldofs = zeros(Int, n) - _celldofs = celldofs!(eldofs, dh, cellnum) - counter = 1 - - for node in cell.nodes - for d in 1:field_dim - data[d, node] = u[_celldofs[counter + offset]] - @debug println(" exporting $(u[_celldofs[counter + offset]]) for dof#$(_celldofs[counter + offset])") - counter += 1 - end - if field_dim == 2 - # paraview requires 3D-data so pad with zero - data[3, node] = 0 - end - end - end - end + data = reshape_to_nodes(dh, u, name) vtk_point_data(vtkfile, data, string(name, suffix)) end return vtkfile -end +end \ No newline at end of file diff --git a/src/FEValues/common_values.jl b/src/FEValues/common_values.jl index bcc0362d94..5e0cf587e0 100644 --- a/src/FEValues/common_values.jl +++ b/src/FEValues/common_values.jl @@ -2,7 +2,7 @@ using Base: @propagate_inbounds -const ScalarValues{dim,T,shape} = Union{CellScalarValues{dim,T,shape},FaceScalarValues{dim,T,shape}} +const ScalarValues{dim,T,shape} = Union{CellScalarValues{dim,T,shape},FaceScalarValues{dim,T,shape},PointScalarValues{dim,T,shape}} const VectorValues{dim,T,shape} = Union{CellVectorValues{dim,T,shape},FaceVectorValues{dim,T,shape}} getnbasefunctions(cv::Values) = size(cv.N, 1) @@ -135,18 +135,19 @@ Base.@pure _valuetype(::VectorValues{dim}, ::AbstractVector{T}) where {dim,T} = # Base.@pure _valuetype(::VectorValues{dim}, ::AbstractVector{Vec{dim,T}}) where {dim,T} = Vec{dim,T} """ - function_scalar_gradient(fe_v::Values{dim}, q_point::Int, u::AbstractVector) + function_gradient(fe_v::Values{dim}, q_point::Int, u::AbstractVector) Compute the gradient of the function in a quadrature point. `u` is a vector with values for the degrees of freedom. For a scalar valued function, `u` contains scalars. For a vector valued function, `u` can be a vector of scalars (for use of `VectorValues`) or `u` can be a vector of `Vec`s (for use with ScalarValues). -The gradient of a scalar function is computed as -``\\mathbf{\\nabla} u(\\mathbf{x}) = \\sum\\limits_{i = 1}^n \\mathbf{\\nabla} N_i (\\mathbf{x}) u_i`` +The gradient of a scalar function or a vector valued function with use of `VectorValues` is computed as +``\\mathbf{\\nabla} u(\\mathbf{x}) = \\sum\\limits_{i = 1}^n \\mathbf{\\nabla} N_i (\\mathbf{x}) u_i`` or +``\\mathbf{\\nabla} u(\\mathbf{x}) = \\sum\\limits_{i = 1}^n \\mathbf{\\nabla} \\mathbf{N}_i (\\mathbf{x}) u_i`` respectively, where ``u_i`` are the nodal values of the function. -For a vector valued function the gradient is computed as -``\\mathbf{\\nabla} \\mathbf{u}(\\mathbf{x}) = \\sum\\limits_{i = 1}^n \\mathbf{\\nabla} N_i (\\mathbf{x}) \\otimes \\mathbf{u}_i`` +For a vector valued function with use of `ScalarValues` the gradient is computed as +``\\mathbf{\\nabla} \\mathbf{u}(\\mathbf{x}) = \\sum\\limits_{i = 1}^n \\mathbf{u}_i \\otimes \\mathbf{\\nabla} N_i (\\mathbf{x})`` where ``\\mathbf{u}_i`` are the nodal values of ``\\mathbf{u}``. """ function function_gradient(fe_v::Values{dim}, q_point::Int, u::AbstractVector{T}, dof_range = eachindex(u)) where {dim,T} diff --git a/src/Ferrite.jl b/src/Ferrite.jl index d8054f97f1..efacec445f 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -5,8 +5,8 @@ using Reexport using LinearAlgebra using SparseArrays - using Base: @propagate_inbounds +using NearestNeighbors include("exports.jl") @@ -43,6 +43,7 @@ include("Quadrature/quadrature.jl") # FEValues include("FEValues/cell_values.jl") include("FEValues/face_values.jl") +include("PointEval/point_values.jl") include("FEValues/common_values.jl") include("FEValues/face_integrals.jl") @@ -67,6 +68,9 @@ include("L2_projection.jl") # Export include("Export/VTK.jl") +# Point Evaluation +include("PointEval/PointEvalHandler.jl") + # Other include("deprecations.jl") diff --git a/src/L2_projection.jl b/src/L2_projection.jl index 68c533e3e2..34725a8be1 100644 --- a/src/L2_projection.jl +++ b/src/L2_projection.jl @@ -85,7 +85,7 @@ function L2Projector( _, vertex_dict, _, _ = __close!(dh) M = _assemble_L2_matrix(fe_values_mass, set, dh) # the "mass" matrix - M_cholesky = cholesky(M) # TODO maybe have a lazy eval instead of precomputing? / JB + M_cholesky = cholesky(M) # For deprecated API fe_values = qr_rhs === nothing ? nothing : @@ -153,7 +153,7 @@ end """ - project(proj::L2Projector, vals::Vector{Vector{T}}}, qr_rhs::QuadratureRule; project_to_nodes=true) + project(proj::L2Projector, vals, qr_rhs::QuadratureRule; project_to_nodes=true) Makes a L2 projection of data `vals` to the nodes of the grid using the projector `proj` (see [`L2Projector`](@ref)). @@ -167,7 +167,9 @@ where ``f`` is the data to project, i.e. `vals`. The data `vals` should be a vector, with length corresponding to number of elements, of vectors, with length corresponding to number of quadrature points per element, matching the number of points in `qr_rhs`. -Example scalar input data: +Alternatively, `vals` can be a matrix, with number of columns corresponding to number of elements, +and number of rows corresponding to number of points in `qr_rhs`. +Example (scalar) input data: ```julia vals = [ [0.44, 0.98, 0.32], # data for quadrature point 1, 2, 3 of element 1 @@ -175,6 +177,14 @@ vals = [ # ... ] ``` +or equivalent in matrix form: +```julia +vals = [ + 0.44 0.29 # ... + 0.98 0.48 # ... + 0.32 0.55 # ... +] +``` Supported data types to project are `Number`s and `AbstractTensor`s. If the parameter `project_to_nodes` is `true`, then the projection returns the values in the order of the mesh nodes @@ -182,7 +192,7 @@ If the parameter `project_to_nodes` is `true`, then the projection returns the v field over the domain, which is useful if one wants to interpolate the projected values. """ function project(proj::L2Projector, - vars::Vector{Vector{T}}, + vars::AbstractVector{<:AbstractVector{T}}, qr_rhs::Union{QuadratureRule,Nothing}=nothing; project_to_nodes::Bool=true) where T <: Union{Number, AbstractTensor} @@ -192,26 +202,29 @@ function project(proj::L2Projector, CellScalarValues(qr_rhs, proj.func_ip, proj.geom_ip) M = T <: AbstractTensor ? length(vars[1][1].data) : 1 - make_T(vals) = T <: AbstractTensor ? T(Tuple(vals)) : vals[1] - projected_vals = _project(vars, proj, fe_values, M) + projected_vals = _project(vars, proj, fe_values, M, T)::Vector{T} if project_to_nodes # NOTE we may have more projected values than verticies in the mesh => not all values are returned nnodes = getnnodes(proj.dh.grid) - reordered_vals = fill(NaN, nnodes, size(projected_vals, 2)) + reordered_vals = fill(convert(T, NaN * zero(T)), nnodes) for node = 1:nnodes - if haskey(proj.node2dof_map, node) - reordered_vals[node, :] = projected_vals[proj.node2dof_map[node], :] + if (k = get(proj.node2dof_map, node, nothing); k !== nothing) + @assert length(k) == 1 + reordered_vals[node] = projected_vals[k[1]] end end - return T[make_T(reordered_vals[i,:]) for i=1:size(reordered_vals, 1)] + return reordered_vals else - # convert back to the original tensor type - return T[make_T(reordered_vals[i,:]) for i=1:size(reordered_vals, 1)] + return projected_vals end end +function project(p::L2Projector, vars::AbstractMatrix, qr_rhs::QuadratureRule; project_to_nodes=true) + # TODO: Random access into vars is required for now, hence the collect + return project(p, collect(eachcol(vars)), qr_rhs; project_to_nodes=project_to_nodes) +end -function _project(vars, proj::L2Projector, fe_values::Values, M::Integer) +function _project(vars, proj::L2Projector, fe_values::Values, M::Integer, ::Type{T}) where {T} # Assemble the multi-column rhs, f = ∭( v ⋅ x̂ )dΩ # The number of columns corresponds to the length of the data-tuple in the tensor x̂. @@ -251,6 +264,9 @@ function _project(vars, proj::L2Projector, fe_values::Values, M::Integer) end # solve for the projected nodal values - return proj.M_cholesky \ f + projected_vals = proj.M_cholesky \ f + # Recast to original input type + make_T(vals) = T <: AbstractTensor ? T(Tuple(vals)) : vals[1] + return T[make_T(x) for x in eachrow(projected_vals)] end diff --git a/src/PointEval/PointEvalHandler.jl b/src/PointEval/PointEvalHandler.jl new file mode 100644 index 0000000000..a6d1a6a3ac --- /dev/null +++ b/src/PointEval/PointEvalHandler.jl @@ -0,0 +1,275 @@ +""" +PointEvalHandler(grid::Grid, points::AbstractVector{Vec{dim,T}}) where {dim, T} + +A `PointEvalHandler` computes the corresponding cell for each point in `points` and its local coordinate within the cell. +""" +PointEvalHandler + +struct PointEvalHandler{G,dim,T<:Real} + grid::G + cells::Vector{Union{Nothing, Int}} + local_coords::Vector{Union{Nothing, Vec{dim,T}}} +end + +function Base.show(io::IO, ::MIME"text/plain", ph::PointEvalHandler) + println(io, typeof(ph)) + println(io, " number of points: ", length(ph.local_coords)) + n_missing = sum(x -> x === nothing, ph.cells) + if n_missing == 0 + print(io, " Found corresponding cell for all points.") + else + print(io, " Could not find corresponding cell for ", n_missing, " points.") + end +end + +function PointEvalHandler(grid::AbstractGrid, points::AbstractVector{Vec{dim,T}}) where {dim, T} + node_cell_dicts = _get_node_cell_map(grid) + cells, local_coords = _get_cellcoords(points, grid, node_cell_dicts) + return PointEvalHandler(grid, cells, local_coords) +end + +function _get_cellcoords(points::AbstractVector{Vec{dim,T}}, grid::AbstractGrid, node_cell_dicts::Dict{C,Dict{Int, Vector{Int}}}) where {dim, T<:Real, C} + + # set up tree structure for finding nearest nodes to points + kdtree = KDTree(reinterpret(Vec{dim,T}, getnodes(grid))) + nearest_nodes, _ = knn(kdtree, points, 3, true) #TODO 3 is a random value, it shouldn't matter because likely the nearest node is the one we want + + cells = Vector{Union{Nothing, Int}}(nothing, length(points)) + local_coords = Vector{Union{Nothing, Vec{dim, T}}}(nothing, length(points)) + + for point_idx in 1:length(points) + cell_found = false + for (idx, (CT, node_cell_dict)) in enumerate(node_cell_dicts) + geom_interpol = default_interpolation(CT) + # loop over points + for node in nearest_nodes[point_idx] + possible_cells = get(node_cell_dict, node, nothing) + possible_cells === nothing && continue # if node is not part of the fieldhandler, try the next node + for cell in possible_cells + cell_coords = getcoordinates(grid, cell) + is_in_cell, local_coord = point_in_cell(geom_interpol, cell_coords, points[point_idx]) + if is_in_cell + cell_found = true + cells[point_idx] = cell + local_coords[point_idx] = local_coord + break + end + end + cell_found && break + end + cell_found && break + end + if !cell_found + @warn("No cell found for point number $point_idx, coordinate: $(points[point_idx]).") + end + end + return cells, local_coords +end + +# check if point is inside a cell based on physical coordinate +function point_in_cell(geom_interpol::Interpolation{dim,shape,order}, cell_coordinates, global_coordinate) where {dim, shape, order} + converged, x_local = find_local_coordinate(geom_interpol, cell_coordinates, global_coordinate) + if converged + return _check_isoparametric_boundaries(shape, x_local), x_local + else + return false, x_local + end +end + +# check if point is inside a cell based on isoparametric coordinate +function _check_isoparametric_boundaries(::Type{RefCube}, x_local::Vec{dim, T}) where {dim, T} + tol = sqrt(eps(T)) + # All in the range [-1, 1] + return all(x -> abs(x) - 1 < tol, x_local) +end + +# check if point is inside a cell based on isoparametric coordinate +function _check_isoparametric_boundaries(::Type{RefTetrahedron}, x_local::Vec{dim, T}) where {dim, T} + tol = sqrt(eps(T)) + # Positive and below the plane 1 - ξx - ξy - ξz + return all(x -> x > -tol, x_local) && sum(x_local) - 1 < tol +end + +# TODO: should we make iteration params optional keyword arguments? +function find_local_coordinate(interpolation, cell_coordinates, global_coordinate) + """ + currently copied verbatim from https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/2 + other than to make J dim x dim rather than 2x2 + """ + dim = length(global_coordinate) + local_guess = zero(Vec{dim}) + n_basefuncs = getnbasefunctions(interpolation) + max_iters = 10 + tol_norm = 1e-10 + converged = false + for iter in 1:max_iters + N = Ferrite.value(interpolation, local_guess) + + global_guess = zero(Vec{dim}) + for j in 1:n_basefuncs + global_guess += N[j] * cell_coordinates[j] + end + residual = global_guess - global_coordinate + if norm(residual) <= tol_norm + converged = true + break + end + dNdξ = Ferrite.derivative(interpolation, local_guess) + J = zero(Tensor{2, dim}) + for j in 1:n_basefuncs + J += cell_coordinates[j] ⊗ dNdξ[j] + end + local_guess -= inv(J) ⋅ residual + end + return converged, local_guess +end + +# return a Dict with a key for each node that contains a vector with the adjacent cells as value +function _get_node_cell_map(grid::AbstractGrid) + cells = getcells(grid) + C = eltype(cells) # possibly abstract + cell_dicts = Dict{Type{<:C}, Dict{Int, Vector{Int}}}() + ctypes = Set{Type{<:C}}(typeof(c) for c in cells) + for ctype in ctypes + cell_dict = cell_dicts[ctype] = Dict{Int,Vector{Int}}() + for (cellidx, cell) in enumerate(cells) + cell isa ctype || continue + for node in cell.nodes + v = get!(Vector{Int}, cell_dict, node) + push!(v, cellidx) + end + end + end + return cell_dicts +end + +""" + get_point_values(ph::PointEvalHandler, proj::L2Projector, dof_values::Vector{T}) where T + get_point_values(ph::PointEvalHandler, dh::AbstractDofHandler, dof_values::Vector{T}, [fieldname::Symbol]) where T + +Return a `Vector{T}` (for a 1-dimensional field) or a `Vector{Vec{fielddim, T}}` (for a vector field) with the field values of field `fieldname` in the points of the `PointEvalHandler`. The `fieldname` can be omitted if only one field is stored in `dh`. The field values are computed based on the `dof_values` and interpolated to the local coordinates by the function interpolation of the corresponding `field` stored in the `AbstractDofHandler` or the `L2Projector`. + +""" +get_point_values + +function get_point_values(ph::PointEvalHandler, proj::L2Projector, dof_vals::AbstractVector) + get_point_values(ph, proj.dh, dof_vals) +end + +function get_point_values(ph::PointEvalHandler, dh::AbstractDofHandler, dof_vals::AbstractVector{T}, + fname::Symbol=find_single_field(dh)) where {T} + fdim = getfielddim(dh, fname) + npoints = length(ph.cells) + # for a scalar field return a Vector of Scalars, for a vector field return a Vector of Vecs + if fdim == 1 + out_vals = fill!(Vector{T}(undef, npoints), NaN * zero(T)) + else + nanv = convert(Vec{fdim,T}, NaN * zero(Vec{fdim,T})) + out_vals = fill!(Vector{Vec{fdim, T}}(undef, npoints), nanv) + end + func_interpolations = get_func_interpolations(dh, fname) + get_point_values!(out_vals, ph, dh, dof_vals, fname, func_interpolations) + return out_vals +end +function find_single_field(dh) + ns = getfieldnames(dh) + if length(ns) != 1 + throw(ArgumentError("multiple fields in DoF handler, must specify which")) + end + return ns[1] +end + +# values in dof-order. They must be obtained from the same DofHandler that was used for constructing the PointEvalHandler +function get_point_values!(out_vals::Vector{T2}, + ph::PointEvalHandler, + dh::MixedDofHandler, + dof_vals::Vector{T}, + fname::Symbol, + func_interpolations + ) where {T2, T} + + # TODO: I don't think this is correct?? + length(dof_vals) == ndofs(dh) || error("You must supply values for all $(ndofs(dh)) dofs.") + + fdim = getfielddim(dh, fname) + + for fh_idx in eachindex(dh.fieldhandlers) + ip = func_interpolations[fh_idx] + if ip !== nothing + dofrange = dof_range(dh.fieldhandlers[fh_idx], fname) + cellset = dh.fieldhandlers[fh_idx].cellset + _get_point_values!(out_vals, dof_vals, ph, dh, ip, cellset, Val(fdim), dofrange) + end + end + return out_vals +end + +function get_point_values!(out_vals::Vector{T2}, + ph::PointEvalHandler, + dh::DofHandler, + dof_vals::Vector{T}, + fname::Symbol, + func_interpolations + ) where {T2, T} + + # TODO: I don't think this is correct?? + length(dof_vals) == ndofs(dh) || error("You must supply values for all $(ndofs(dh)) dofs.") + + fdim = getfielddim(dh, fname) + dofrange = dof_range(dh, fname) + _get_point_values!(out_vals, dof_vals, ph, dh, func_interpolations[1], nothing, Val(fdim), dofrange) + return out_vals +end + +# function barrier with concrete type of interpolation +function _get_point_values!( + out_vals::Vector{T2}, + dof_vals::Vector{T}, + ph::PointEvalHandler, + dh::AbstractDofHandler, + ip::Interpolation, + cellset::Union{Nothing, Set{Int}}, + fdim::Val{fielddim}, + dofrange::AbstractRange{Int}, + ) where {T2,T,fielddim} + + # extract variables + local_coords = ph.local_coords + # preallocate some stuff specific to this cellset + pv = PointScalarValues(first(local_coords), ip) + first_cell = cellset === nothing ? 1 : first(cellset) + cell_dofs = Vector{Int}(undef, ndofs_per_cell(dh, first_cell)) + + # compute point values + for pointid in eachindex(ph.cells) + cellid = ph.cells[pointid] + cellid ===nothing && continue # next point if no cell was found for this one + cellset !== nothing && (cellid ∈ cellset || continue) # no need to check the cellset for a regular DofHandler + celldofs!(cell_dofs, dh, ph.cells[pointid]) + reinit!(pv, local_coords[pointid], ip) + @inbounds @views dof_vals_reshaped = _change_format(fdim, dof_vals[cell_dofs[dofrange]]) + out_vals[pointid] = function_value(pv, 1, dof_vals_reshaped) + end + return out_vals +end + +################################################################################################################### +# utils + +# reshape dof_values based on fielddim +_change_format(::Val{1}, dof_values::AbstractVector{T}) where T = dof_values +_change_format(::Val{fielddim}, dof_values::AbstractVector{T}) where {fielddim, T} = reinterpret(Vec{fielddim, T}, dof_values) + +get_func_interpolations(dh::DH, fieldname) where DH<:DofHandler = [getfieldinterpolation(dh, find_field(dh, fieldname))] +function get_func_interpolations(dh::DH, fieldname) where DH<:MixedDofHandler + func_interpolations = Union{Interpolation,Nothing}[] + for fh in dh.fieldhandlers + j = findfirst(i -> i === fieldname, getfieldnames(fh)) + if j === nothing + push!(func_interpolations, missing) + else + push!(func_interpolations, fh.fields[j].interpolation) + end + end + return func_interpolations +end diff --git a/src/PointEval/point_values.jl b/src/PointEval/point_values.jl new file mode 100644 index 0000000000..5d5f081b85 --- /dev/null +++ b/src/PointEval/point_values.jl @@ -0,0 +1,47 @@ + +struct PointScalarValues{dim,T<:Real,refshape<:AbstractRefShape} <: CellValues{dim,T,refshape} + N::Vector{T} +end + +function Base.show(io::IO, ::MIME"text/plain", pv::PointScalarValues) + print(io, "$(typeof(pv)) with $(getnbasefunctions(pv)) shape functions.") +end + +function PointScalarValues(quad_rule::QuadratureRule, func_interpol::Interpolation) + PointScalarValues(Float64, quad_rule, func_interpol) +end + +function PointScalarValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_interpol::Interpolation) where {dim,T,shape<:AbstractRefShape} + + length(getweights(quad_rule)) == 1 || error("PointScalarValues supports only a single point.") + + # Function interpolation + n_func_basefuncs = getnbasefunctions(func_interpol) + N = fill(zero(T) * T(NaN), n_func_basefuncs) + + ξ = quad_rule.points[1] + for i in 1:n_func_basefuncs + N[i] = value(func_interpol, i, ξ) + end + + PointScalarValues{dim,T,shape}(N) +end + +# PointScalarValues only have one quadrature point anyways +function PointScalarValues(coord::Vec{dim,T}, ip::Interpolation{dim, refshape}) where {dim,refshape,T} + qr = QuadratureRule{dim,refshape,T}([one(T)], [coord]) + return PointScalarValues(qr, ip) +end + +# allow to use function_value with any +Base.@pure _valuetype(::PointScalarValues{dim}, ::Vector{T}) where {dim, T<:AbstractTensor} = T + +# allow on-the-fly updating +function reinit!(pv::PointScalarValues{dim,T,refshape}, coord::Vec{dim,T}, func_interpol::Interpolation{dim,refshape,order}) where {dim,T,refshape,order} + n_func_basefuncs = getnbasefunctions(func_interpol) + for i in 1:n_func_basefuncs + pv.N[i] = value(func_interpol, i, coord) + end + return pv +end +# TODO: need a show method for PointScalarValues \ No newline at end of file diff --git a/src/exports.jl b/src/exports.jl index b70cd7437b..7f7d82ce11 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -100,6 +100,7 @@ export MixedDofHandler, FieldHandler, Field, + reshape_to_nodes, # Constraints ConstraintHandler, @@ -132,4 +133,8 @@ export # L2 Projection project, - L2Projector + L2Projector, + +# Point Evaluation + PointEvalHandler, + get_point_values diff --git a/test/checksums2.sha1 b/test/checksums2.sha1 new file mode 100644 index 0000000000..f68fcec74f --- /dev/null +++ b/test/checksums2.sha1 @@ -0,0 +1,2 @@ +9e5923d04f15fd2cc3abf06ce1725d251c1f09f6 +2d11788599fb980fe0279e32c5e7720d6e1087ba diff --git a/test/runtests.jl b/test/runtests.jl index 1662eccdb6..c5f8f13ffe 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,7 @@ include("test_constraints.jl") include("test_grid_dofhandler_vtk.jl") include("test_mixeddofhandler.jl") include("test_l2_projection.jl") +include("test_pointevaluation.jl") # include("test_notebooks.jl") include("test_apply_rhs.jl") include("test_examples.jl") diff --git a/test/test_dofs.jl b/test/test_dofs.jl index 9e832b299b..eb28ae31b4 100644 --- a/test/test_dofs.jl +++ b/test/test_dofs.jl @@ -91,4 +91,18 @@ close!(dh) @test celldofs(dh,1) == collect(1:18) @test celldofs(dh,2) == [2, 19, 20, 3, 21, 22, 23, 6, 24, 11, 25, 26, 12, 27, 28, 29, 15, 30] +# test reshape_to_nodes +## DofHandler +mesh = generate_grid(Quadrilateral, (1,1)) +dh = DofHandler(mesh) +push!(dh, :v, 2) +push!(dh, :s, 1) +close!(dh) + +u = [1.1, 1.2, 2.1, 2.2, 4.1, 4.2, 3.1, 3.2, 1.3, 2.3, 4.3, 3.3] + +s_nodes = reshape_to_nodes(dh, u, :s) +@test s_nodes ≈ [i+0.3 for i=1:4]' +v_nodes = reshape_to_nodes(dh, u, :v) +@test v_nodes ≈ [i==3 ? 0.0 : j+i/10 for i=1:3, j=1:4] end \ No newline at end of file diff --git a/test/test_examples.jl b/test/test_examples.jl index ed3bd17cbb..fa702de9d9 100644 --- a/test/test_examples.jl +++ b/test/test_examples.jl @@ -30,4 +30,12 @@ module TestQuasiIncompressibleHyperElasticity include(joinpath(@__DIR__, "../docs/src/literate/quasi_incompressible_hyperelasticity.jl")) end end -end \ No newline at end of file +end + +# module TestNavierStokesDiffeqIntegration +# mktempdir() do dir +# cd(dir) do +# include(joinpath(@__DIR__, "../docs/src/literate/ns_vs_diffeq.jl")) +# end +# end +# end \ No newline at end of file diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 65fe207746..f53dab0247 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/test/test_grid_dofhandler_vtk.jl @@ -95,9 +95,59 @@ end # of testset close(csio) +@testset "vtk tensor export" begin + # open files + checksums_file_tensors = joinpath(dirname(@__FILE__), "checksums2.sha1") + if OVERWRITE_CHECKSUMS + csio = open(checksums_file_tensors, "w") + else + csio = open(checksums_file_tensors, "r") + end + + # 3D grid + grid = generate_grid(Hexahedron, (1,1,1)) + + sym_tensor_data = [SymmetricTensor{2,3}(ntuple(i->i, 6)) for j=1.0:8.0] + tensor_data = [Tensor{2,3}(ntuple(i->i, 9)) for j=1.0:8.0] + vector_data = [Vec{3}(ntuple(i->i, 3)) for j=1:8] + + filename_3d = "test_vtk_3d" + vtk_grid(filename_3d, grid) do vtk_file + vtk_point_data(vtk_file, sym_tensor_data, "symmetric tensor") + vtk_point_data(vtk_file, tensor_data, "tensor") + vtk_point_data(vtk_file, vector_data, "vector") + end + + # 2D grid + grid = generate_grid(Quadrilateral, (1,1)) + + sym_tensor_data = [SymmetricTensor{2,2}(ntuple(i->i, 3)) for j=1.0:4.0] + tensor_data = [Tensor{2,2}(ntuple(i->i, 4)) for j=1.0:4.0] + tensor_data_1D = [SymmetricTensor{2,1}(ntuple(i->i, 1)) for j=1.0:4.0] + vector_data = [Vec{2}(ntuple(i->i, 2)) for j=1:4] + + filename_2d = "test_vtk_2d" + vtk_grid(filename_2d, grid) do vtk_file + vtk_point_data(vtk_file, sym_tensor_data, "symmetric tensor") + vtk_point_data(vtk_file, tensor_data, "tensor") + vtk_point_data(vtk_file, tensor_data_1D, "tensor_1d") + vtk_point_data(vtk_file, vector_data, "vector") + end -# right = Vec{2, Float64}(ntuple(x->1.5, dim)) -# left = -right + # test the shas of the files + files = [filename_3d, filename_2d] + for filename in files + sha = bytes2hex(open(SHA.sha1, filename*".vtu")) + if OVERWRITE_CHECKSUMS + write(csio, sha, "\n") + else + @test chomp(readline(csio)) == sha + rm(filename*".vtu") + end + end + + close(csio) +end @testset "Grid utils" begin diff --git a/test/test_l2_projection.jl b/test/test_l2_projection.jl index bbbbe27f54..6121955c10 100644 --- a/test/test_l2_projection.jl +++ b/test/test_l2_projection.jl @@ -39,14 +39,16 @@ function test_projection(order, refshape) # Now recover the nodal values using a L2 projection. proj = L2Projector(ip, grid; geom_ip=ip_geom) point_vars = project(proj, qp_values, qr) + qp_values_matrix = reduce(hcat, qp_values) + point_vars_2 = project(proj, qp_values_matrix, qr) ## Old API with fe values as first arg proj2 = @test_deprecated L2Projector(cv, ip, grid) - point_vars_2 = @test_deprecated project(qp_values, proj2) + point_vars_3 = @test_deprecated project(qp_values, proj2) ## Old API with qr as first arg proj3 = @test_deprecated L2Projector(qr, ip, grid) - point_vars_3 = @test_deprecated project(qp_values, proj3) + point_vars_4 = @test_deprecated project(qp_values, proj3) - @test point_vars ≈ point_vars_2 ≈ point_vars_3 + @test point_vars ≈ point_vars_2 ≈ point_vars_3 ≈ point_vars_4 if order == 1 # A linear approximation can not recover a quadratic solution, @@ -72,7 +74,12 @@ function test_projection(order, refshape) # Tensor f_tensor(x) = Tensor{2,2,Float64}((f(x),2*f(x),3*f(x),4*f(x))) qp_values = analytical(f_tensor) + qp_values_matrix = reduce(hcat, qp_values)::Matrix point_vars = project(proj, qp_values, qr) + point_vars_2 = project(proj, qp_values_matrix, qr) + + @test point_vars ≈ point_vars_2 + if order == 1 ae = [Tensor{2,2,Float64}((f_approx(i),2*f_approx(i),3*f_approx(i),4*f_approx(i))) for i in 1:4] elseif order == 2 @@ -83,7 +90,12 @@ function test_projection(order, refshape) # SymmetricTensor f_stensor(x) = SymmetricTensor{2,2,Float64}((f(x),2*f(x),3*f(x))) qp_values = analytical(f_stensor) + qp_values_matrix = reduce(hcat, qp_values) point_vars = project(proj, qp_values, qr) + point_vars_2 = project(proj, qp_values_matrix, qr) + + @test point_vars ≈ point_vars_2 + if order == 1 ae = [SymmetricTensor{2,2,Float64}((f_approx(i),2*f_approx(i),3*f_approx(i))) for i in 1:4] elseif order == 2 @@ -132,6 +144,7 @@ function test_projection_mixedgrid() ae = compute_vertex_values(mesh, f) # analytical values qp_values = [[f(spatial_coordinate(cv, qp, xe)) for qp in 1:getnquadpoints(cv)]] + qp_values_matrix = reduce(hcat, qp_values) # Now recover the nodal values using a L2 projection. # Assume f would only exist on the first cell, we project it to the nodes of the @@ -139,16 +152,18 @@ function test_projection_mixedgrid() # nodes that do not belong to the 1st cell proj = L2Projector(ip, mesh; geom_ip=ip_geom, set=1:1) point_vars = project(proj, qp_values, qr) + point_vars_2 = project(proj, qp_values_matrix, qr) ## Old API with fe values as first arg proj = @test_deprecated L2Projector(cv, ip, mesh, 1:1) - point_vars_2 = @test_deprecated project(qp_values, proj) + point_vars_3 = @test_deprecated project(qp_values, proj) ## Old API with qr as first arg proj = @test_deprecated L2Projector(qr, ip, mesh, 1:1) - point_vars_3 = @test_deprecated project(qp_values, proj) + point_vars_4 = @test_deprecated project(qp_values, proj) # In the nodes of the 1st cell we should recover the field for node in mesh.cells[1].nodes - @test ae[node] ≈ point_vars[node] ≈ point_vars_2[node] ≈ point_vars_3[node] + @test ae[node] ≈ point_vars[node] ≈ point_vars_2[node] ≈ point_vars_3[node] ≈ + point_vars_4[node] end # in all other nodes we should have NaNs @@ -157,14 +172,37 @@ function test_projection_mixedgrid() @test isnan(point_vars[node][d1, d2]) @test isnan(point_vars_2[node][d1, d2]) @test isnan(point_vars_3[node][d1, d2]) + @test isnan(point_vars_4[node][d1, d2]) end end end +function test_node_reordering() + grid = generate_grid(Quadrilateral, (1, 1), Vec((0.,0.)), Vec((2.,2.))) + dim = 2 + ip = Lagrange{dim, RefCube, 2}() + ip_geo = Lagrange{dim, RefCube,1}() + qr = QuadratureRule{dim, RefCube}(3) + cv = CellScalarValues(qr, ip, ip_geo) + + f(x) = x[1]+x[2] + + qp_values = [[f(spatial_coordinate(cv, qp, getcoordinates(cell))) for qp in 1:getnquadpoints(cv)] for cell in CellIterator(grid)] + + projector = L2Projector(ip, grid) + projected_vals_nodes = project(projector, qp_values, qr) + projected_vals_dofs = project(projector, qp_values, qr; project_to_nodes=false) + + tol = 1e-12 + @test all(projected_vals_nodes - [0.0, 2.0, 2.0, 4.0] .< tol) + @test all(projected_vals_dofs - [0., 2., 4., 2., 1., 3., 3., 1., 2.] .< tol) +end + @testset "Test L2-Projection" begin test_projection(1, RefCube) test_projection(1, RefTetrahedron) test_projection(2, RefCube) test_projection(2, RefTetrahedron) test_projection_mixedgrid() + test_node_reordering() end diff --git a/test/test_mixeddofhandler.jl b/test/test_mixeddofhandler.jl index 73d0e5f544..a45b43ea1f 100644 --- a/test/test_mixeddofhandler.jl +++ b/test/test_mixeddofhandler.jl @@ -434,6 +434,55 @@ function test_field_on_subdomain() @test_throws ErrorException Ferrite.find_field(dh.fieldhandlers[1], :s) end +function test_reshape_to_nodes() + + # 5_______6 + # |\ | + # | \ | + # 3______\4 + # | | + # | | + # 1_______2 + + nodes = [Node((0.0, 0.0)), + Node((1.0, 0.0)), + Node((0.0, 1.0)), + Node((1.0, 1.0)), + Node((0.0, 2.0)), + Node((1.0, 2.0))] + cells = Ferrite.AbstractCell[Quadrilateral((1,2,4,3)), + Triangle((3,4,6)), + Triangle((3,6,5))] + mesh = Grid(cells, nodes) + addcellset!(mesh, "quads", Set{Int}((1,))) + addcellset!(mesh, "tris", Set{Int}((2, 3))) + + ip_quad = Lagrange{2,RefCube,1}() + ip_tri = Lagrange{2,RefTetrahedron,1}() + + dh = MixedDofHandler(mesh) + field_v_tri = Field(:v, ip_tri, 2) # vector field :v everywhere + fh_tri = FieldHandler([field_v_tri], getcellset(mesh, "tris")) + push!(dh, fh_tri) + field_v_quad = Field(:v, ip_quad, 2) + field_s_quad = Field(:s, ip_quad, 1) # scalar field :s only on quad + fh_quad = FieldHandler([field_v_quad, field_s_quad], getcellset(mesh, "quads")) + push!(dh, fh_quad) + close!(dh) + + u = collect(1.:16.) + + s_nodes = reshape_to_nodes(dh, u, :s) + @test s_nodes[1:4] ≈ [13., 14., 16., 15.] + @test all(isnan.(s_nodes[5:6])) + v_nodes = reshape_to_nodes(dh, u, :v) + @test v_nodes ≈ hcat( [9., 10., 0.], + [11., 12., 0.], + [1., 2., 0.], + [3., 4., 0.], + [7., 8., 0.], + [5., 6., 0.]) +end function test_mixed_grid_show() grid = get_2d_grid() @@ -548,5 +597,7 @@ end test_mixed_grid_show(); test_subparametric_quad(); test_subparametric_triangle(); + test_reshape_to_nodes() + test_mixed_grid_show() test_separate_fields_on_separate_domains(); end diff --git a/test/test_pointevaluation.jl b/test/test_pointevaluation.jl new file mode 100644 index 0000000000..be0543fa15 --- /dev/null +++ b/test/test_pointevaluation.jl @@ -0,0 +1,257 @@ +function scalar_field() + # isoparametric approximation + mesh = generate_grid(QuadraticQuadrilateral, (20, 20)) + f(x) = x[1]^2 + nodal_vals = [f(p.x) for p in mesh.nodes] + + ip_f = Lagrange{2,RefCube,2}() # function interpolation + ip_g = Lagrange{2,RefCube,2}() # geometry interpolation + + # compute values in quadrature points + qr = QuadratureRule{2, RefCube}(3) # exactly approximate quadratic field + cv = CellScalarValues(qr, ip_f, ip_g) + qp_vals = [Vector{Float64}(undef, getnquadpoints(cv)) for _ in 1:getncells(mesh)] + for cellid in eachindex(mesh.cells) + xe = getcoordinates(mesh, cellid) + reinit!(cv, xe) + for qp in 1:getnquadpoints(cv) + qp_vals[cellid][qp] = f(spatial_coordinate(cv, qp, xe)) + end + end + + # do a L2Projection for getting values in dofs + projector = L2Projector(ip_f, mesh) + projector_vals = project(projector, qp_vals, qr; project_to_nodes=false) + + # points where we want to retrieve field values + points = [Vec((x, 0.52)) for x in range(0.0; stop=1.0, length=100)] + + # set up PointEvalHandler and retrieve values + ph = PointEvalHandler(mesh, points) + vals = get_point_values(ph, projector, projector_vals) + @test f.(points) ≈ vals + + # alternatively retrieve vals from nodal values TODO: make this work? + # vals = get_point_values(ph, nodal_vals) + # @test f.(points) ≈ vals +end + +function vector_field() + ## vector field + # isoparametric approximation + mesh = generate_grid(QuadraticQuadrilateral, (20, 20)) + f(x) = Vec((x[1]^2, x[1])) + nodal_vals = [f(p.x) for p in mesh.nodes] + + ip_f = Lagrange{2,RefCube,2}() # function interpolation + ip_g = Lagrange{2,RefCube,2}() # geometry interpolation + + # compute values in quadrature points + qr = QuadratureRule{2, RefCube}(3) # exactly approximate quadratic field + cv = CellScalarValues(qr, ip_f, ip_g) + qp_vals = [Vector{Vec{2,Float64}}(undef, getnquadpoints(cv)) for i=1:getncells(mesh)] + for cellid in eachindex(mesh.cells) + xe = getcoordinates(mesh, cellid) + reinit!(cv, xe) + for qp in 1:getnquadpoints(cv) + qp_vals[cellid][qp] = f(spatial_coordinate(cv, qp, xe)) + end + end + + # do a L2Projection for getting values in dofs + projector = L2Projector(ip_f, mesh) + projector_vals = project(projector, qp_vals, qr; project_to_nodes=false) + # TODO: project_to_nodes should probably return dof values and not Vecs for vector fields + # projector_vals = convert(Vector{Float64}, reinterpret(Float64, projector_vals)) + + # points where we want to retrieve field values + points = [Vec((x, 0.52)) for x in range(0.0; stop=1.0, length=100)] + + # set up PointEvalHandler and retrieve values + ph = PointEvalHandler(mesh, points) + vals = get_point_values(ph, projector, projector_vals) + @test f.(points) ≈ vals + + # alternatively retrieve vals from nodal values# TODO + # vals = get_point_values(ph, nodal_vals) + # @test f.(points) ≈ vals +end + +function superparametric() + # superparametric approximation + mesh = generate_grid(Quadrilateral, (20, 20)) + f(x) = x*x[1] + ip_f = Lagrange{2,RefCube,2}() # function interpolation + ip_g = Lagrange{2,RefCube,1}() # geometry interpolation + + # compute values in quadrature points + qr = QuadratureRule{2, RefCube}(3) # exactly approximate quadratic field + cv = CellScalarValues(qr, ip_f, ip_g) + qp_vals = [Vector{Vec{2,Float64}}(undef, getnquadpoints(cv)) for i=1:getncells(mesh)] + for cellid in eachindex(mesh.cells) + xe = getcoordinates(mesh, cellid) + reinit!(cv, xe) + for qp in 1:getnquadpoints(cv) + qp_vals[cellid][qp] = f(spatial_coordinate(cv, qp, xe)) + end + end + + # do a L2Projection for getting values in dofs + projector = L2Projector(ip_f, mesh) + projector_vals = project(projector, qp_vals, qr; project_to_nodes=false) + + # points where we want to retrieve field values + points = [Vec((x, 0.52)) for x in range(0.0; stop=1.0, length=100)] + + # set up PointEvalHandler and retrieve values + ph = PointEvalHandler(mesh, points) + vals = get_point_values(ph, projector, projector_vals) + + # can recover a quadratic field by a quadratic approximation + @test f.(points) ≈ vals +end + +function dofhandler() + mesh = generate_grid(Quadrilateral, (2,2)) + dof_vals = [1., 2., 5., 4., 3., 6., 8., 7., 9.] + points = [node.x for node in mesh.nodes] # same as nodes + + dh = DofHandler(mesh) + push!(dh, :s, 1) # a scalar field + close!(dh) + + ph = PointEvalHandler(mesh, points) + vals = get_point_values(ph, dh, dof_vals, :s) + @test vals ≈ 1.0:9.0 + + # TODO + # vals = get_point_values(ph, collect(1.0:9.0)) + # @test vals ≈ 1.0:9.0 +end + +function mixed_grid() + ## Mixed grid where not all cells have the same fields + + # 5_______6 + # |\ | + # | \ | + # 3______\4 + # | | + # | | + # 1_______2 + + nodes = [Node((0.0, 0.0)), + Node((1.0, 0.0)), + Node((0.0, 1.0)), + Node((1.0, 1.0)), + Node((0.0, 2.0)), + Node((1.0, 2.0))] + + cells = Ferrite.AbstractCell[Quadrilateral((1,2,4,3)), + Triangle((3,4,6)), + Triangle((3,6,5))] + + mesh = Grid(cells, nodes) + addcellset!(mesh, "quads", Set{Int}((1,))) + addcellset!(mesh, "tris", Set{Int}((2, 3))) + + ip_quad = Lagrange{2,RefCube,1}() + ip_tri = Lagrange{2,RefTetrahedron,1}() + + f(x) = x[1] + + # compute values in quadrature points for quad + qr = QuadratureRule{2, RefCube}(2) + cv = CellScalarValues(qr, ip_quad) + qp_vals_quads = [Vector{Float64}(undef, getnquadpoints(cv)) for cell in getcellset(mesh, "quads")] + for (local_cellid, global_cellid) in enumerate(getcellset(mesh, "quads")) + xe = getcoordinates(mesh, global_cellid) + reinit!(cv, xe) + for qp in 1:getnquadpoints(cv) + qp_vals_quads[local_cellid][qp] = f(spatial_coordinate(cv, qp, xe)) + end + end + + # construct projector + projector = L2Projector(ip_quad, mesh; set=getcellset(mesh, "quads")) + + points = [Vec((x, 2x)) for x in range(0.0; stop=1.0, length=10)] + + # first alternative: L2Projection to dofs + projector_values = project(projector, qp_vals_quads, qr; project_to_nodes = false) + ph = PointEvalHandler(mesh, points) + vals = get_point_values(ph, projector, projector_values) + @test vals[1:5] ≈ f.(points[1:5]) + @test all(isnan, vals[6:end]) + # TODO + # # second alternative: L2Projection to nodes + # nodal_vals = project(projector, qp_vals_quads, qr; project_to_nodes = true) + # vals = get_point_values(ph, nodal_vals) + # @test vals[1:5] ≈ f.(points[1:5]) + # @test all(isnan.(vals[6:end])) + + + # second alternative: assume a vector field :v + dh = MixedDofHandler(mesh) + field = Field(:v, ip_quad, 2) + fh_quad = FieldHandler([field], getcellset(mesh, "quads")) + push!(dh, fh_quad) + field = Field(:v, ip_tri, 2) + fh_tri = FieldHandler([field], getcellset(mesh, "tris")) + push!(dh, fh_tri) + close!(dh) + + dof_vals = [1., 1., 2., 2., 4., 4., 3., 3., 6., 6., 5., 5.] + points = [node.x for node in mesh.nodes] + ph = PointEvalHandler(mesh, points) + vals = get_point_values(ph, dh, dof_vals, :v) + @test vals == [Vec((i, i)) for i=1.0:6.0] +end + +function oneD() + # isoparametric approximation + mesh = generate_grid(Line, (2,)) + f(x) = x[1] + nodal_vals = [f(p.x) for p in mesh.nodes] + + ip_f = Lagrange{1,RefCube,1}() # function interpolation + ip_g = Lagrange{1,RefCube,1}() # geometry interpolation + + # compute values in quadrature points + qr = QuadratureRule{1, RefCube}(2) + cv = CellScalarValues(qr, ip_f, ip_g) + qp_vals = [Vector{Float64}(undef, getnquadpoints(cv)) for i=1:getncells(mesh)] + for cellid in eachindex(mesh.cells) + xe = getcoordinates(mesh, cellid) + reinit!(cv, xe) + for qp in 1:getnquadpoints(cv) + qp_vals[cellid][qp] = f(spatial_coordinate(cv, qp, xe)) + end + end + + # do a L2Projection for getting values in dofs + projector = L2Projector(ip_f, mesh) + projector_values = project(projector, qp_vals, qr; project_to_nodes=false) + + # points where we want to retrieve field values + points = [Vec((x,)) for x in range(-1.0; stop=1.0, length=5)] + + # set up PointEvalHandler and retrieve values + ph = PointEvalHandler(mesh, points) + vals = get_point_values(ph, projector, projector_values) + @test f.(points) ≈ vals + + # alternatively retrieve vals from nodal values + # TODO + # vals = get_point_values(ph, nodal_vals) + # @test f.(points) ≈ vals +end + +@testset "PointEvalHandler" begin + scalar_field() + vector_field() + dofhandler() + superparametric() + mixed_grid() + oneD() +end