From 3fb52077279531ceb65c3bf2463dd122c4e0d074 Mon Sep 17 00:00:00 2001 From: David Rollin Date: Mon, 6 Jul 2026 13:54:43 +0200 Subject: [PATCH 01/12] Fix quadratic ip on linear geo --- .github/workflows/CI.yml | 1 + Manifest.toml | 357 +++++++++++++++++++++++++++++++++ Project.toml | 4 +- src/cellvalues.jl | 7 +- test/test_cellvalues.jl | 2 + test/test_example_cohesion.jl | 16 +- test/test_example_diffusion.jl | 16 +- test/test_examples.jl | 56 +++--- 8 files changed, 412 insertions(+), 47 deletions(-) create mode 100644 Manifest.toml diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 5cf205e..e132fbd 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -24,6 +24,7 @@ jobs: matrix: version: - '1.10' + - '1.12' - 'nightly' os: - ubuntu-latest diff --git a/Manifest.toml b/Manifest.toml new file mode 100644 index 0000000..0d038e2 --- /dev/null +++ b/Manifest.toml @@ -0,0 +1,357 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.12.5" +manifest_format = "2.0" +project_hash = "6c5f3321cfce62ff245e91360f1dd30d1bbae77d" + +[[deps.AbstractTrees]] +git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177" +uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" +version = "0.4.5" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +version = "1.11.0" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +version = "1.11.0" + +[[deps.CodecZlib]] +deps = ["TranscodingStreams", "Zlib_jll"] +git-tree-sha1 = "962834c22b66e32aa10f7611c08c8ca4e20749a9" +uuid = "944b1d66-785c-5afd-91f1-9de20f533193" +version = "0.7.8" + +[[deps.CommonSubexpressions]] +deps = ["MacroTools"] +git-tree-sha1 = "cda2cfaebb4be89c9084adaca7dd7333369715c5" +uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" +version = "0.3.1" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.3.0+1" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +version = "1.11.0" + +[[deps.DiffResults]] +deps = ["StaticArraysCore"] +git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621" +uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +version = "1.1.0" + +[[deps.DiffRules]] +deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] +git-tree-sha1 = "79a2aca180a85c690c58a020d47b426954b590f8" +uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" +version = "1.16.0" + +[[deps.Distances]] +deps = ["LinearAlgebra", "Statistics", "StatsAPI"] +git-tree-sha1 = "c7e3a542b999843086e2f29dac96a618c105be1d" +uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +version = "0.10.12" + + [deps.Distances.extensions] + DistancesChainRulesCoreExt = "ChainRulesCore" + DistancesSparseArraysExt = "SparseArrays" + + [deps.Distances.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.DocStringExtensions]] +git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.9.5" + +[[deps.EnumX]] +git-tree-sha1 = "c49898e8438c828577f04b92fc9368c388ac783c" +uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56" +version = "1.0.7" + +[[deps.Ferrite]] +deps = ["EnumX", "ForwardDiff", "LinearAlgebra", "NearestNeighbors", "OrderedCollections", "Preferences", "Reexport", "SparseArrays", "Tensors", "WriteVTK"] +git-tree-sha1 = "f61701b7c7b7a5feee1c8a82b390c63dfa4c9da6" +uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992" +version = "1.4.1" + + [deps.Ferrite.extensions] + FerriteBlockArrays = "BlockArrays" + FerriteMetis = "Metis" + FerriteSparseMatrixCSR = "SparseMatricesCSR" + + [deps.Ferrite.weakdeps] + BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" + Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" + SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" + +[[deps.FerriteInterfaceElements]] +deps = ["Ferrite", "OrderedCollections"] +path = "." +uuid = "76cf2c14-24ca-44b2-b86c-6046d3cfcf82" +version = "1.0.0" + +[[deps.FillArrays]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "2f979084d1e13948a3352cf64a25df6bd3b4dca3" +uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" +version = "1.16.0" + + [deps.FillArrays.extensions] + FillArraysPDMatsExt = "PDMats" + FillArraysSparseArraysExt = "SparseArrays" + FillArraysStaticArraysExt = "StaticArrays" + FillArraysStatisticsExt = "Statistics" + + [deps.FillArrays.weakdeps] + PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[deps.ForwardDiff]] +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"] +git-tree-sha1 = "2c5d0b0e12088cde2cf84afb2784415b1ea3dfee" +uuid = "f6369f11-7733-5829-9624-2563aa707210" +version = "1.4.1" +weakdeps = ["StaticArrays"] + + [deps.ForwardDiff.extensions] + ForwardDiffStaticArraysExt = "StaticArrays" + +[[deps.IrrationalConstants]] +git-tree-sha1 = "b2d91fe939cae05960e760110b328288867b5758" +uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" +version = "0.2.6" + +[[deps.JLLWrappers]] +deps = ["Artifacts", "Preferences"] +git-tree-sha1 = "7204148362dafe5fe6a273f855b8ccbe4df8173e" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.8.0" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +version = "1.11.0" + +[[deps.Libiconv_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "be484f5c92fad0bd8acfef35fe017900b0b73809" +uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" +version = "1.18.0+0" + +[[deps.LightXML]] +deps = ["Libdl", "XML2_jll"] +git-tree-sha1 = "aa971a09f0f1fe92fe772713a564aa48abe510df" +uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179" +version = "0.9.3" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +version = "1.12.0" + +[[deps.LogExpFunctions]] +deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "bba2d9aa057d8f126415de240573e86a8f39d2a1" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "1.0.1" + + [deps.LogExpFunctions.extensions] + LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" + LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" + LogExpFunctionsInverseFunctionsExt = "InverseFunctions" + + [deps.LogExpFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + +[[deps.MacroTools]] +git-tree-sha1 = "1e0228a030642014fe5cfe68c2c0a818f9e3f522" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.16" + +[[deps.NaNMath]] +deps = ["OpenLibm_jll"] +git-tree-sha1 = "dbd2e8cd2c1c27f0b584f6661b4309609c5a685e" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "1.1.4" + +[[deps.NearestNeighbors]] +deps = ["AbstractTrees", "Distances", "StaticArrays"] +git-tree-sha1 = "ca562494d657e2b69191e440e9c28b9692d67944" +uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce" +version = "0.4.28" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.29+0" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.7+0" + +[[deps.OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1346c9208249809840c91b26703912dff463d335" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.6+0" + +[[deps.OrderedCollections]] +git-tree-sha1 = "94ba93778373a53bfd5a0caaf7d809c445292ff4" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.8.2" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "edbeefc7a4889f528644251bdb5fc9ab5348bc2c" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.3.4" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "8b770b60760d4451834fe79dd483e318eee709c4" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.5.2" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +version = "1.11.0" + +[[deps.Random]] +deps = ["SHA"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +version = "1.11.0" + +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "1.2.2" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.SIMD]] +deps = ["PrecompileTools"] +git-tree-sha1 = "e24dc23107d426a096d3eae6c165b921e74c18e4" +uuid = "fdea26ae-647d-5447-a871-4b548cad5224" +version = "3.7.2" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" +version = "1.11.0" + +[[deps.SparseArrays]] +deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +version = "1.12.0" + +[[deps.SpecialFunctions]] +deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "6547cbdd8ce32efba0d21c5a40fa96d1a3548f9f" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "2.8.0" + + [deps.SpecialFunctions.extensions] + SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" + + [deps.SpecialFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] +git-tree-sha1 = "246a8bb2e6667f832eea063c3a56aef96429a3db" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.9.18" + + [deps.StaticArrays.extensions] + StaticArraysChainRulesCoreExt = "ChainRulesCore" + StaticArraysStatisticsExt = "Statistics" + + [deps.StaticArrays.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "6ab403037779dae8c514bad259f32a447262455a" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.4" + +[[deps.Statistics]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "ae3bb1eb3bba077cd276bc5cfc337cc65c3075c0" +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +version = "1.11.1" +weakdeps = ["SparseArrays"] + + [deps.Statistics.extensions] + SparseArraysExt = ["SparseArrays"] + +[[deps.StatsAPI]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "178ed29fd5b2a2cfc3bd31c13375ae925623ff36" +uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" +version = "1.8.0" + +[[deps.SuiteSparse_jll]] +deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] +uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" +version = "7.8.3+2" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.Tensors]] +deps = ["ForwardDiff", "LinearAlgebra", "PrecompileTools", "SIMD", "StaticArrays", "Statistics"] +git-tree-sha1 = "28b845db2855a43f0228dfb2e1c69cac08e53553" +uuid = "48a634ad-e948-5137-8d70-aa71f2a747f4" +version = "1.17.1" + +[[deps.TranscodingStreams]] +git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742" +uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" +version = "0.11.3" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" +version = "1.11.0" + +[[deps.VTKBase]] +git-tree-sha1 = "c2d0db3ef09f1942d08ea455a9e252594be5f3b6" +uuid = "4004b06d-e244-455f-a6ce-a5f9919cc534" +version = "1.0.1" + +[[deps.WriteVTK]] +deps = ["Base64", "CodecZlib", "FillArrays", "LightXML", "TranscodingStreams", "VTKBase"] +git-tree-sha1 = "073f2ae23cc1aa11510772d4c156435cdb8d7087" +uuid = "64499a7a-5c06-52f2-abe2-ccb03c286192" +version = "1.22.0" + +[[deps.XML2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"] +git-tree-sha1 = "3f3315d89fc954a28f5b471bce698ed6e27481be" +uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" +version = "2.15.3+0" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.3.1+2" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.15.0+0" diff --git a/Project.toml b/Project.toml index 219b9ae..318ff96 100644 --- a/Project.toml +++ b/Project.toml @@ -8,9 +8,9 @@ Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" [compat] -Ferrite = "1.0" +Ferrite = "1.4" OrderedCollections = "1" -julia = "1.10" +julia = "1.12" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/cellvalues.jl b/src/cellvalues.jl index 98d1409..9b37e2a 100755 --- a/src/cellvalues.jl +++ b/src/cellvalues.jl @@ -105,14 +105,15 @@ Ferrite.shape_value_type(cv::InterfaceCellValues) = shape_value_type(cv.here) Ferrite.shape_gradient_type(cv::InterfaceCellValues) = shape_gradient_type(cv.here) -Ferrite.reinit!(cv::InterfaceCellValues, cc::CellCache) = reinit!(cv, cc.coords) # TODO: Needed? +Ferrite.reinit!(cv::InterfaceCellValues, cc::CellCache) = reinit!(cv, cc.coords) function Ferrite.reinit!(cv::InterfaceCellValues{CV}, x::AbstractVector{Vec{sdim,T}}) where {sdim, T, CV} - x_here = @view x[cv.base_indices_here] + n_coords_per_side = length(x) ÷ 2 + x_here = @view x[cv.base_indices_here[1:n_coords_per_side]] reinit!(cv.here, x_here) if ! (cv.here === cv.there) - x_there = @view x[cv.base_indices_there] + x_there = @view x[cv.base_indices_there[1:n_coords_per_side]] reinit!(cv.there, x_there) end diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 3186fd4..fee05f0 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -35,6 +35,8 @@ @test all(abs.(function_gradient_average(cv, qp, u)) .≤ 1e-14) @test all(abs.(function_gradient_jump(cv, qp, u)) .≤ 1e-14) @test getdetJdV_average(cv, qp) == (getdetJdV(cv.here, qp) + getdetJdV(cv.there, qp)) / 2 + n = @allocated function_value_jump(cv, qp, u) + @test n == 0 end qr = QuadratureRule{RefTriangle}(2) diff --git a/test/test_example_cohesion.jl b/test/test_example_cohesion.jl index 16c948a..f2a745c 100644 --- a/test/test_example_cohesion.jl +++ b/test/test_example_cohesion.jl @@ -31,14 +31,16 @@ function assemble_test_element_cohesion!(Kₑ::Matrix, cv::InterfaceCellValues, end end -@testset "In $(dim)D for order $(order)" for (order, dim, get_grid, get_ip) in ( - (1, 2, prepare_interface_test_grid_2D, prepare_interface_test_vector_interpolation_2D), - (2, 2, prepare_interface_test_grid_2D, prepare_interface_test_vector_interpolation_2D), - (1, 3, prepare_interface_test_grid_3D, prepare_interface_test_vector_interpolation_3D), - (2, 3, prepare_interface_test_grid_3D, prepare_interface_test_vector_interpolation_3D) +@testset "In $(dim)D for order (ip=$(iporder), grid=($gridorder))" for (iporder, gridorder, dim, get_grid, get_ip) in ( + (1, 1, 2, prepare_interface_test_grid_2D, prepare_interface_test_vector_interpolation_2D), + (2, 2, 2, prepare_interface_test_grid_2D, prepare_interface_test_vector_interpolation_2D), + (2, 1, 2, prepare_interface_test_grid_2D, prepare_interface_test_vector_interpolation_2D), + (1, 1, 3, prepare_interface_test_grid_3D, prepare_interface_test_vector_interpolation_3D), + (2, 2, 3, prepare_interface_test_grid_3D, prepare_interface_test_vector_interpolation_3D), + (2, 1, 3, prepare_interface_test_grid_3D, prepare_interface_test_vector_interpolation_3D) ) - grid = get_grid(order) - ip, cv = get_ip(order) + grid = get_grid(gridorder) + ip, cv = get_ip(iporder, gridorder) dh = DofHandler(grid) add!(SubDofHandler(dh, Set{Int}([1])), :u, ip.left) diff --git a/test/test_example_diffusion.jl b/test/test_example_diffusion.jl index ded48d3..023154d 100644 --- a/test/test_example_diffusion.jl +++ b/test/test_example_diffusion.jl @@ -28,14 +28,16 @@ function assemble_test_element_diffusion!(Kₑ::Matrix, cv::InterfaceCellValues) end end -@testset "In $(dim)D for order $(order)" for (order, dim, get_grid, get_ip) in ( - (1, 2, prepare_interface_test_grid_2D, prepare_interface_test_scalar_interpolation_2D), - (2, 2, prepare_interface_test_grid_2D, prepare_interface_test_scalar_interpolation_2D), - (1, 3, prepare_interface_test_grid_3D, prepare_interface_test_scalar_interpolation_3D), - (2, 3, prepare_interface_test_grid_3D, prepare_interface_test_scalar_interpolation_3D) +@testset "In $(dim)D for order (ip=$(iporder), grid=$(gridorder))" for (iporder, gridorder, dim, get_grid, get_ip) in ( + (1, 1, 2, prepare_interface_test_grid_2D, prepare_interface_test_scalar_interpolation_2D), + (2, 2, 2, prepare_interface_test_grid_2D, prepare_interface_test_scalar_interpolation_2D), + (2, 1, 2, prepare_interface_test_grid_2D, prepare_interface_test_scalar_interpolation_2D), + (1, 1, 3, prepare_interface_test_grid_3D, prepare_interface_test_scalar_interpolation_3D), + (2, 2, 3, prepare_interface_test_grid_3D, prepare_interface_test_scalar_interpolation_3D), + (2, 1, 3, prepare_interface_test_grid_3D, prepare_interface_test_scalar_interpolation_3D) ) - grid = get_grid(order) - ip, cv = get_ip(order) + grid = get_grid(gridorder) + ip, cv = get_ip(iporder, gridorder) dh = DofHandler(grid) add!(SubDofHandler(dh, Set{Int}([1])), :c, ip.left) diff --git a/test/test_examples.jl b/test/test_examples.jl index 8fa1778..1253cd7 100644 --- a/test/test_examples.jl +++ b/test/test_examples.jl @@ -62,47 +62,47 @@ function prepare_interface_test_grid_3D(order::Integer) return grid end -function prepare_interface_test_scalar_interpolation_2D(order::Integer) - ip = ( left = Lagrange{RefQuadrilateral, order}(), - interface = InterfaceCellInterpolation(Lagrange{RefLine, order}()), - right = Lagrange{RefQuadrilateral, order}()) +function prepare_interface_test_scalar_interpolation_2D(iporder::Integer, gridorder::Int) + ip = ( left = Lagrange{RefQuadrilateral, iporder}(), + interface = InterfaceCellInterpolation(Lagrange{RefLine, iporder}()), + right = Lagrange{RefQuadrilateral, iporder}()) - cv = ( left = CellValues(QuadratureRule{RefQuadrilateral}(4), ip.left, ip.left), - interface = InterfaceCellValues(QuadratureRule{RefLine}(4), ip.interface), - right = CellValues(QuadratureRule{RefQuadrilateral}(4), ip.right, ip.right)) + cv = ( left = CellValues(QuadratureRule{RefQuadrilateral}(4), ip.left, Lagrange{RefQuadrilateral, gridorder}()), + interface = InterfaceCellValues(QuadratureRule{RefLine}(4), ip.interface, InterfaceCellInterpolation(Lagrange{RefLine, gridorder}())^2), + right = CellValues(QuadratureRule{RefQuadrilateral}(4), ip.right, Lagrange{RefQuadrilateral, gridorder}())) return ip, cv end -function prepare_interface_test_scalar_interpolation_3D(order::Integer) - ip = ( left = Lagrange{RefHexahedron, order}(), - interface = InterfaceCellInterpolation(Lagrange{RefQuadrilateral, order}()), - right = Lagrange{RefHexahedron, order}()) +function prepare_interface_test_scalar_interpolation_3D(iporder::Integer, gridorder::Int) + ip = ( left = Lagrange{RefHexahedron, iporder}(), + interface = InterfaceCellInterpolation(Lagrange{RefQuadrilateral, iporder}()), + right = Lagrange{RefHexahedron, iporder}()) - cv = ( left = CellValues(QuadratureRule{RefHexahedron}(4), ip.left, ip.left), - interface = InterfaceCellValues(QuadratureRule{RefQuadrilateral}(4), ip.interface), - right = CellValues(QuadratureRule{RefHexahedron}(4), ip.right, ip.right)) + cv = ( left = CellValues(QuadratureRule{RefHexahedron}(4), ip.left, Lagrange{RefHexahedron, gridorder}()), + interface = InterfaceCellValues(QuadratureRule{RefQuadrilateral}(4), ip.interface, InterfaceCellInterpolation(Lagrange{RefQuadrilateral, gridorder}())^3), + right = CellValues(QuadratureRule{RefHexahedron}(4), ip.right, Lagrange{RefHexahedron, gridorder}())) return ip, cv end -function prepare_interface_test_vector_interpolation_2D(order::Integer) - ip = ( left = Lagrange{RefQuadrilateral, order}()^2, - interface = InterfaceCellInterpolation(Lagrange{RefLine, order}())^2, - right = Lagrange{RefQuadrilateral, order}()^2) +function prepare_interface_test_vector_interpolation_2D(iporder::Integer, gridorder::Int) + ip = ( left = Lagrange{RefQuadrilateral, iporder}()^2, + interface = InterfaceCellInterpolation(Lagrange{RefLine, iporder}())^2, + right = Lagrange{RefQuadrilateral, iporder}()^2) - cv = ( left = CellValues(QuadratureRule{RefQuadrilateral}(4), ip.left, ip.left), - interface = InterfaceCellValues(QuadratureRule{RefLine}(4), ip.interface), - right = CellValues(QuadratureRule{RefQuadrilateral}(4), ip.right, ip.right)) + cv = ( left = CellValues(QuadratureRule{RefQuadrilateral}(4), ip.left, Lagrange{RefQuadrilateral, gridorder}()^2), + interface = InterfaceCellValues(QuadratureRule{RefLine}(4), ip.interface, InterfaceCellInterpolation(Lagrange{RefLine, gridorder}())^2), + right = CellValues(QuadratureRule{RefQuadrilateral}(4), ip.right, Lagrange{RefQuadrilateral, gridorder}()^2)) return ip, cv end -function prepare_interface_test_vector_interpolation_3D(order::Integer) - ip = ( left = Lagrange{RefHexahedron, order}()^3, - interface = InterfaceCellInterpolation(Lagrange{RefQuadrilateral, order}())^3, - right = Lagrange{RefHexahedron, order}()^3) +function prepare_interface_test_vector_interpolation_3D(iporder::Integer, gridorder::Int) + ip = ( left = Lagrange{RefHexahedron, iporder}()^3, + interface = InterfaceCellInterpolation(Lagrange{RefQuadrilateral, iporder}())^3, + right = Lagrange{RefHexahedron, iporder}()^3) - cv = ( left = CellValues(QuadratureRule{RefHexahedron}(4), ip.left, ip.left), - interface = InterfaceCellValues(QuadratureRule{RefQuadrilateral}(4), ip.interface), - right = CellValues(QuadratureRule{RefHexahedron}(4), ip.right, ip.right)) + cv = ( left = CellValues(QuadratureRule{RefHexahedron}(4), ip.left, Lagrange{RefHexahedron, gridorder}()^3), + interface = InterfaceCellValues(QuadratureRule{RefQuadrilateral}(4), ip.interface, InterfaceCellInterpolation(Lagrange{RefQuadrilateral, gridorder}())^3), + right = CellValues(QuadratureRule{RefHexahedron}(4), ip.right, Lagrange{RefHexahedron, gridorder}()^3)) return ip, cv end From ad93dabbabe21f437b917a5658b886f6d43dde08 Mon Sep 17 00:00:00 2001 From: David Rollin Date: Mon, 6 Jul 2026 14:21:13 +0200 Subject: [PATCH 02/12] Update grid tests to OrderedSet --- test/runtests.jl | 3 ++- test/test_grid.jl | 32 ++++++++++++++++---------------- 2 files changed, 18 insertions(+), 17 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 706d241..9da8c54 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,6 @@ using FerriteInterfaceElements -using FerriteInterfaceElements.Ferrite +using Ferrite +using OrderedCollections using Test @testset "FerriteInterfaceElements.jl" begin diff --git a/test/test_grid.jl b/test/test_grid.jl index 0485a70..f401098 100644 --- a/test/test_grid.jl +++ b/test/test_grid.jl @@ -3,8 +3,8 @@ # |\ c6 |\ c8 | |\ c6 | c|\ c8 | # | \ | \ | | \ |11| \ | # | c5 \| c7 \| | c5 \| | c7 \| - # 4 ___ 5 ___ 6 ---> 13___12__11___10 - # |\ c2 |\ c4 | | c10 \/ c9 | + # 4 ___ 5 ___ 6 ---> 11___10__13___12 + # |\ c2 |\ c4 | | c9 \/ c10 | # | \ | \ | 4 ___ 5 ___ 6 # | c1 \| c3 \| |\ c2 | \ c4 | # 1 ___ 2 ___ 3 | \ | \ | @@ -12,14 +12,14 @@ # 1 ____ 2 ____ 3 # grid = generate_grid(Triangle, (2,2)) - addcellset!(grid, "bottom", Set((1,2,3,4))) - addcellset!(grid, "topleft", Set((5,6))) - addcellset!(grid, "topright", Set((7,8))) + addcellset!(grid, "bottom", OrderedSet((1,2,3,4))) + addcellset!(grid, "topleft", OrderedSet((5,6))) + addcellset!(grid, "topright", OrderedSet((7,8))) domain_names = ["bottom", "topleft", "topright"] new_grid = insert_interfaces(grid, domain_names) - for (nodeid, duplicate_nodeid) in [(4,13), (5,12), (5,11), (11,12), (6,10), (8,14)] + for (nodeid, duplicate_nodeid) in [(4,11), (5,10), (5,13), (6,12), (8,14), (10,13)] @test new_grid.nodes[nodeid] == new_grid.nodes[duplicate_nodeid] end @@ -28,21 +28,21 @@ Triangle((2, 5, 4)), Triangle((2, 3, 5)), Triangle((3, 6, 5)), - Triangle((13, 12, 7)), - Triangle((12, 8, 7)), - Triangle((11, 10, 14)), - Triangle((10, 9, 14)), - InterfaceCell(Line((6, 5)), Line((10, 11))), - InterfaceCell(Line((5, 4)), Line((12, 13))), - InterfaceCell(Line((12, 8)), Line((11, 14))) + Triangle((11, 10, 7)), + Triangle((10, 8, 7)), + Triangle((13, 12, 14)), + Triangle((12, 9, 14)), + InterfaceCell(Line((5, 4)), Line((10, 11))), + InterfaceCell(Line((6, 5)), Line((12, 13))), + InterfaceCell(Line((10, 8)), Line((13, 14))) ] end @testset "Inserting interfaces in 3D" begin grid = generate_grid(Hexahedron, (2,2,1)) - addcellset!(grid, "bottomleft", Set((1,))) - addcellset!(grid, "topleft", Set((3,))) - addcellset!(grid, "right", Set((2,4))) + addcellset!(grid, "bottomleft", OrderedSet((1,))) + addcellset!(grid, "topleft", OrderedSet((3,))) + addcellset!(grid, "right", OrderedSet((2,4))) domain_names = ["bottomleft", "topleft", "right"] new_grid = insert_interfaces(grid, domain_names) From 003843aa9da8d6ba0b1bdc1550ef173d4b94583d Mon Sep 17 00:00:00 2001 From: David Rollin Date: Mon, 6 Jul 2026 14:51:29 +0200 Subject: [PATCH 03/12] Make grid tests independent of cell order --- test/test_grid.jl | 50 ++++++++++++++++++++++++----------------------- 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/test/test_grid.jl b/test/test_grid.jl index f401098..061444c 100644 --- a/test/test_grid.jl +++ b/test/test_grid.jl @@ -22,20 +22,21 @@ for (nodeid, duplicate_nodeid) in [(4,11), (5,10), (5,13), (6,12), (8,14), (10,13)] @test new_grid.nodes[nodeid] == new_grid.nodes[duplicate_nodeid] end - - @test new_grid.cells == [ - Triangle((1, 2, 4)), - Triangle((2, 5, 4)), - Triangle((2, 3, 5)), - Triangle((3, 6, 5)), - Triangle((11, 10, 7)), - Triangle((10, 8, 7)), - Triangle((13, 12, 14)), - Triangle((12, 9, 14)), - InterfaceCell(Line((5, 4)), Line((10, 11))), - InterfaceCell(Line((6, 5)), Line((12, 13))), - InterfaceCell(Line((10, 8)), Line((13, 14))) - ] + for cell in [Triangle((1, 2, 4)), + Triangle((2, 5, 4)), + Triangle((2, 3, 5)), + Triangle((3, 6, 5)), + Triangle((11, 10, 7)), + Triangle((10, 8, 7)), + Triangle((13, 12, 14)), + Triangle((12, 9, 14))] + @test cell in new_grid.cells[1:8] + end + for cell in [InterfaceCell(Line((5, 4)), Line((10, 11))), + InterfaceCell(Line((6, 5)), Line((12, 13))), + InterfaceCell(Line((10, 8)), Line((13, 14)))] + @test cell in new_grid.cells[9:11] + end end @testset "Inserting interfaces in 3D" begin @@ -50,16 +51,17 @@ end for (nodeid, duplicate_nodeid) in [(13, 25), (14,26), (14,21), (21,26), (11,22), (17,27), (4,24), (5,23), (5,20), (20,23), (8, 28), (2,19)] @test new_grid.nodes[nodeid] == new_grid.nodes[duplicate_nodeid] end - - @test new_grid.cells == [ - Hexahedron((1, 2, 5, 4, 10, 11, 14, 13)), - Hexahedron((19, 3, 6, 20, 22, 12, 15, 21)), - Hexahedron((24, 23, 28, 7, 25, 26, 27, 16)), - Hexahedron((20, 6, 9, 8, 21, 15, 18, 17)), - InterfaceCell(Quadrilateral((2, 5, 14, 11)), Quadrilateral((19, 20, 21, 22))), - InterfaceCell(Quadrilateral((5, 4, 13, 14)), Quadrilateral((23, 24, 25, 26))), - InterfaceCell(Quadrilateral((20, 21, 17, 8)), Quadrilateral((23, 26, 27, 28))) - ] + for cell in [Hexahedron((1, 2, 5, 4, 10, 11, 14, 13)), + Hexahedron((19, 3, 6, 20, 22, 12, 15, 21)), + Hexahedron((24, 23, 28, 7, 25, 26, 27, 16)), + Hexahedron((20, 6, 9, 8, 21, 15, 18, 17))] + @test cell in new_grid.cells[1:4] + end + for cell in [InterfaceCell(Quadrilateral((2, 5, 14, 11)), Quadrilateral((19, 20, 21, 22))), + InterfaceCell(Quadrilateral((5, 4, 13, 14)), Quadrilateral((23, 24, 25, 26))), + InterfaceCell(Quadrilateral((20, 21, 17, 8)), Quadrilateral((23, 26, 27, 28)))] + @test cell in new_grid.cells[5:7] + end # More complex example with rough test grid = generate_grid(Tetrahedron, (10,10,10), Vec((-0.5,-0.5,-0.5)), Vec((0.5,0.5,0.5))) From a3a8a1be5025ce1c01fd48c628fb8e932e311158 Mon Sep 17 00:00:00 2001 From: David Rollin Date: Mon, 6 Jul 2026 15:05:45 +0200 Subject: [PATCH 04/12] Improve grid test to automatically find node pairs --- test/test_grid.jl | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/test/test_grid.jl b/test/test_grid.jl index 061444c..bcc9ec6 100644 --- a/test/test_grid.jl +++ b/test/test_grid.jl @@ -19,8 +19,16 @@ domain_names = ["bottom", "topleft", "topright"] new_grid = insert_interfaces(grid, domain_names) - for (nodeid, duplicate_nodeid) in [(4,11), (5,10), (5,13), (6,12), (8,14), (10,13)] - @test new_grid.nodes[nodeid] == new_grid.nodes[duplicate_nodeid] + nodepairs = Set{Tuple{Int,Int}}() + for cell in new_grid.cells + if cell isa InterfaceCell + for (n, m) in zip(cell.here.nodes, cell.there.nodes) + push!(nodepairs, (n,m)) + end + end + end + for (n, m) in nodepairs + @test new_grid.nodes[n] == new_grid.nodes[m] end for cell in [Triangle((1, 2, 4)), Triangle((2, 5, 4)), @@ -48,8 +56,16 @@ end domain_names = ["bottomleft", "topleft", "right"] new_grid = insert_interfaces(grid, domain_names) - for (nodeid, duplicate_nodeid) in [(13, 25), (14,26), (14,21), (21,26), (11,22), (17,27), (4,24), (5,23), (5,20), (20,23), (8, 28), (2,19)] - @test new_grid.nodes[nodeid] == new_grid.nodes[duplicate_nodeid] + nodepairs = Set{Tuple{Int,Int}}() + for cell in new_grid.cells + if cell isa InterfaceCell + for (n, m) in zip(cell.here.nodes, cell.there.nodes) + push!(nodepairs, (n,m)) + end + end + end + for (n, m) in nodepairs + @test new_grid.nodes[n] == new_grid.nodes[m] end for cell in [Hexahedron((1, 2, 5, 4, 10, 11, 14, 13)), Hexahedron((19, 3, 6, 20, 22, 12, 15, 21)), From affe8313a46f3154bb76389346464077448c64ed Mon Sep 17 00:00:00 2001 From: David Rollin Date: Mon, 6 Jul 2026 16:48:08 +0200 Subject: [PATCH 05/12] Changed grid tests to be less specific --- test/test_grid.jl | 47 +++++++++++++++++++++++------------------------ 1 file changed, 23 insertions(+), 24 deletions(-) diff --git a/test/test_grid.jl b/test/test_grid.jl index bcc9ec6..82a8c4e 100644 --- a/test/test_grid.jl +++ b/test/test_grid.jl @@ -1,3 +1,22 @@ +function test_grid_data(old::Grid, new::Grid) + # Check if all InterfaceCell's contain new nodes on at least one side and that node pairs have the same coordinates + n = length(old.nodes) + nodepairs = Set{Tuple{Int,Int}}() + for cell in new.cells + if cell isa InterfaceCell + h, t = cell.here.nodes, cell.there.nodes + @test (any( h .≤ n ) && all( t .> n )) || (all( h .> n ) && any( t .≤ n )) + for (n, m) in zip(h, t) + push!(nodepairs, (n,m)) + end + end + end + for (n, m) in nodepairs + @test new.nodes[n] == new.nodes[m] + end + return nothing +end + @testset "Inserting interfaces in 2D" begin # 7 ___ 8 ___ 9 7 ___ 8__14___ 9 # |\ c6 |\ c8 | |\ c6 | c|\ c8 | @@ -18,18 +37,8 @@ domain_names = ["bottom", "topleft", "topright"] new_grid = insert_interfaces(grid, domain_names) - - nodepairs = Set{Tuple{Int,Int}}() - for cell in new_grid.cells - if cell isa InterfaceCell - for (n, m) in zip(cell.here.nodes, cell.there.nodes) - push!(nodepairs, (n,m)) - end - end - end - for (n, m) in nodepairs - @test new_grid.nodes[n] == new_grid.nodes[m] - end + test_grid_data(grid, new_grid) + for cell in [Triangle((1, 2, 4)), Triangle((2, 5, 4)), Triangle((2, 3, 5)), @@ -55,18 +64,8 @@ end domain_names = ["bottomleft", "topleft", "right"] new_grid = insert_interfaces(grid, domain_names) - - nodepairs = Set{Tuple{Int,Int}}() - for cell in new_grid.cells - if cell isa InterfaceCell - for (n, m) in zip(cell.here.nodes, cell.there.nodes) - push!(nodepairs, (n,m)) - end - end - end - for (n, m) in nodepairs - @test new_grid.nodes[n] == new_grid.nodes[m] - end + test_grid_data(grid, new_grid) + for cell in [Hexahedron((1, 2, 5, 4, 10, 11, 14, 13)), Hexahedron((19, 3, 6, 20, 22, 12, 15, 21)), Hexahedron((24, 23, 28, 7, 25, 26, 27, 16)), From 681ec6ae9c1a27ca98eb6075abfc413f35717dbb Mon Sep 17 00:00:00 2001 From: David Rollin Date: Mon, 6 Jul 2026 16:53:12 +0200 Subject: [PATCH 06/12] Remove tests causing issues --- test/test_grid.jl | 28 ---------------------------- 1 file changed, 28 deletions(-) diff --git a/test/test_grid.jl b/test/test_grid.jl index 82a8c4e..f3d6925 100644 --- a/test/test_grid.jl +++ b/test/test_grid.jl @@ -38,22 +38,6 @@ end domain_names = ["bottom", "topleft", "topright"] new_grid = insert_interfaces(grid, domain_names) test_grid_data(grid, new_grid) - - for cell in [Triangle((1, 2, 4)), - Triangle((2, 5, 4)), - Triangle((2, 3, 5)), - Triangle((3, 6, 5)), - Triangle((11, 10, 7)), - Triangle((10, 8, 7)), - Triangle((13, 12, 14)), - Triangle((12, 9, 14))] - @test cell in new_grid.cells[1:8] - end - for cell in [InterfaceCell(Line((5, 4)), Line((10, 11))), - InterfaceCell(Line((6, 5)), Line((12, 13))), - InterfaceCell(Line((10, 8)), Line((13, 14)))] - @test cell in new_grid.cells[9:11] - end end @testset "Inserting interfaces in 3D" begin @@ -65,18 +49,6 @@ end domain_names = ["bottomleft", "topleft", "right"] new_grid = insert_interfaces(grid, domain_names) test_grid_data(grid, new_grid) - - for cell in [Hexahedron((1, 2, 5, 4, 10, 11, 14, 13)), - Hexahedron((19, 3, 6, 20, 22, 12, 15, 21)), - Hexahedron((24, 23, 28, 7, 25, 26, 27, 16)), - Hexahedron((20, 6, 9, 8, 21, 15, 18, 17))] - @test cell in new_grid.cells[1:4] - end - for cell in [InterfaceCell(Quadrilateral((2, 5, 14, 11)), Quadrilateral((19, 20, 21, 22))), - InterfaceCell(Quadrilateral((5, 4, 13, 14)), Quadrilateral((23, 24, 25, 26))), - InterfaceCell(Quadrilateral((20, 21, 17, 8)), Quadrilateral((23, 26, 27, 28)))] - @test cell in new_grid.cells[5:7] - end # More complex example with rough test grid = generate_grid(Tetrahedron, (10,10,10), Vec((-0.5,-0.5,-0.5)), Vec((0.5,0.5,0.5))) From e44cb0c02095cc85b08d95f646da9bae9057ffc2 Mon Sep 17 00:00:00 2001 From: David Rollin Date: Tue, 7 Jul 2026 09:13:12 +0200 Subject: [PATCH 07/12] Improve grid testing --- test/test_grid.jl | 65 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 64 insertions(+), 1 deletion(-) diff --git a/test/test_grid.jl b/test/test_grid.jl index f3d6925..bfcac6f 100644 --- a/test/test_grid.jl +++ b/test/test_grid.jl @@ -2,8 +2,10 @@ function test_grid_data(old::Grid, new::Grid) # Check if all InterfaceCell's contain new nodes on at least one side and that node pairs have the same coordinates n = length(old.nodes) nodepairs = Set{Tuple{Int,Int}}() - for cell in new.cells + interfacecellids = Set{Int}() + for (i, cell) in enumerate(new.cells) if cell isa InterfaceCell + push!(interfacecellids, i) h, t = cell.here.nodes, cell.there.nodes @test (any( h .≤ n ) && all( t .> n )) || (all( h .> n ) && any( t .≤ n )) for (n, m) in zip(h, t) @@ -13,6 +15,67 @@ function test_grid_data(old::Grid, new::Grid) end for (n, m) in nodepairs @test new.nodes[n] == new.nodes[m] + end + # Reconstruct cellsets based on interfaces + nodesets = [Set{Int}(), Set{Int}(), Set{Int}(), Set{Int}()] + cell = new.cells[pop!(interfacecellids)] + for (n, m) in zip(cell.here.nodes, cell.there.nodes) + push!(nodesets[1], n) + push!(nodesets[2], m) + end + # Find nodes on same side of interfaces + while ! isempty(interfacecellids) + cell = new.cells[pop!(interfacecellids)] + for nodes in (cell.here.nodes, cell.there.nodes) + # Search for a set with node overlap + foundconnection = false + for set in nodesets + if any( [n in set for n in nodes] ) + foundconnection = true + for n in nodes + push!(set, n) + end + break + end + end + # If no node overlap, add nodes to empty set + if ! foundconnection + for set in nodesets + if isempty(set) + for n in nodes + push!(set, n) + end + break + end + end + end + # Merge sets with overlap + for i in 1:length(nodesets) + for j in i+1:length(nodesets) + if ! isempty( nodesets[i] ∩ nodesets[j] ) + union!(nodesets[i], nodesets[j]) + empty!(nodesets[j]) + end + end + end + end + end + filter!(s -> ! isempty(s), nodesets) # One set should be empty + # Collect cells according to node sets + cellsets = [OrderedSet{Int}(), OrderedSet{Int}(), OrderedSet{Int}()] + for (cellid, cell) in enumerate(new.cells) + if ! (cell isa InterfaceCell) + for (i, set) in enumerate(nodesets) + if any([n in set for n in cell.nodes]) + push!(cellsets[i], cellid) + break + end + end + end + end + # Test reconstructed cellsets + for set in cellsets + @test any( [sort(set) == sort(oldset) for oldset in values(old.cellsets)] ) end return nothing end From c79c80f35848ec9a479859010b33b34de876f8d4 Mon Sep 17 00:00:00 2001 From: David Rollin Date: Tue, 7 Jul 2026 10:44:45 +0200 Subject: [PATCH 08/12] Fix interface cell insertion for mixed grid (2D) --- src/grid.jl | 37 +++++++++++++++++++++---------------- test/test_grid.jl | 28 ++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+), 16 deletions(-) diff --git a/src/grid.jl b/src/grid.jl index 9af4bd6..ede67b8 100644 --- a/src/grid.jl +++ b/src/grid.jl @@ -1,11 +1,10 @@ ###################################################################### # Inserting cells into a grid ###################################################################### - """ insert_interfaces(grid, domain_names; topology=ExclusiveTopology(grid)) -Return a new grid with `InterfaceCell`s inserted betweenthe domains defined by `domain_names`. +Return a new grid with `InterfaceCell`s inserted between the domains defined by `domain_names`. The new grid provides additional cell sets. The set `"interfaces"` contains all new `InterfaceCell`s and two sets are provided for each combination of domain names: `"domain1-domain2-interface"` and `"domain2-domain1-interface"` both using the same `Set`. @@ -27,7 +26,7 @@ function insert_interfaces(grid, domain_names; topology=ExclusiveTopology(grid)) nodes = copy(grid.nodes) cells_generic = Vector{Ferrite.AbstractCell}(grid.cells) # copies - while !isempty(cellsets) + while ! isempty(cellsets) name, cellset = pop!(cellsets) for cellid in cellset cell = getcells(grid, cellid) @@ -53,7 +52,7 @@ function insert_interfaces(grid, domain_names; topology=ExclusiveTopology(grid)) # generate a new node for (_name, _cellset) push!(nodes, nodes[nodeid]) node_mapping[_name][nodeid] = length(nodes) # main node - elseif isnothing(new_nodeid) && !isnothing(_new_nodeid) + elseif isnothing(new_nodeid) && ! isnothing(_new_nodeid) # node has been duplicated at least once, so it already has a main grain # generate a new node for (name, cellset) push!(nodes, nodes[nodeid]) @@ -63,8 +62,7 @@ function insert_interfaces(grid, domain_names; topology=ExclusiveTopology(grid)) new_nodeids = Tuple(node_mapping[name][i] for i in facetnodeids) _new_nodeids = Tuple(node_mapping[_name][i] for i in facetnodeids) # generate new cell - @assert typeof(cell) == typeof(getcells(grid, _cellid)) - interface_cell = create_interface_cell(typeof(cell), new_nodeids, _new_nodeids) + interface_cell = create_interface_cell(typeof(cell), typeof(getcells(grid, _cellid)), new_nodeids, _new_nodeids) push!(cells_generic, interface_cell) push!(interfacesets["$(name)-$(_name)-interface"], length(cells_generic)) end @@ -94,22 +92,29 @@ function insert_interfaces(grid, domain_names; topology=ExclusiveTopology(grid)) end """ - create_interface_cell(::Type{C}, nodes_here, nodes_there) where {C} + create_interface_cell(::Type{C₁}, ::Type{C₂}, nodes_here, nodes_there) where {C₁,C₂} Return a suitable `InterfaceCell` connecting the facets with `nodes_here` and `nodes_there`. """ -function create_interface_cell(::Type{C}, nodes_here, nodes_there) where {C} - Cbase = get_interface_base_cell_type(C) +function create_interface_cell(::Type{C₁}, ::Type{C₂}, nodes_here, nodes_there) where {C₁,C₂} + Cbase = get_interface_base_cell_type(C₁, C₂) return InterfaceCell(Cbase(nodes_here), Cbase(nodes_there)) end + + """ - get_interface_base_cell_type(::Type{<:AbstractCell}) + get_interface_base_cell_type(::Type{<:AbstractCell}, ::Type{<:AbstractCell}) Return a suitable base type for connecting two cells of given type with an `InterfaceCell`. """ -get_interface_base_cell_type(::Type{Triangle}) = Line -get_interface_base_cell_type(::Type{QuadraticTriangle}) = QuadraticLine -get_interface_base_cell_type(::Type{Quadrilateral}) = Line -get_interface_base_cell_type(::Type{QuadraticQuadrilateral}) = QuadraticLine -get_interface_base_cell_type(::Type{Tetrahedron}) = Triangle -get_interface_base_cell_type(::Type{Hexahedron}) = Quadrilateral +get_interface_base_cell_type(::Type{Triangle}, ::Type{Triangle}) = Line +get_interface_base_cell_type(::Type{QuadraticTriangle}, ::Type{QuadraticTriangle}) = QuadraticLine +get_interface_base_cell_type(::Type{Quadrilateral}, ::Type{Quadrilateral}) = Line +get_interface_base_cell_type(::Type{QuadraticQuadrilateral}, ::Type{QuadraticQuadrilateral}) = QuadraticLine +get_interface_base_cell_type(::Type{Tetrahedron}, ::Type{Tetrahedron}) = Triangle +get_interface_base_cell_type(::Type{Hexahedron}, ::Type{Hexahedron}) = Quadrilateral + +get_interface_base_cell_type(::Type{Triangle}, ::Type{Quadrilateral}) = Line +get_interface_base_cell_type(::Type{Quadrilateral}, ::Type{Triangle}) = Line +get_interface_base_cell_type(::Type{QuadraticTriangle}, ::Type{QuadraticQuadrilateral}) = QuadraticLine +get_interface_base_cell_type(::Type{QuadraticQuadrilateral}, ::Type{QuadraticTriangle}) = QuadraticLine diff --git a/test/test_grid.jl b/test/test_grid.jl index bfcac6f..6015711 100644 --- a/test/test_grid.jl +++ b/test/test_grid.jl @@ -133,3 +133,31 @@ end @test length(getcells(newgrid, "1-3-interface")) == 192 @test length(getcells(newgrid, "2-3-interface")) == 48 end + +@testset "Inserting interfaces in 2D for a mixed grid" begin + # 7 ___ 8 ___ 9 7 ___ 8__14___ 9 + # | | | | | | | + # | c4 | c5 | | c4 |c8| c5 | + # | | | | | | | + # 4 ___ 5 ___ 6 ---> 11___10__13___12 + # |\ c2 | | | c6 \/ c7 | + # | \ | c3 | 4 ___ 5 ___ 6 + # | c1 \| | |\ c2 | | + # 1 ___ 2 ___ 3 | \ | c3 | + # | c1 \ | | + # 1 ____ 2 ____ 3 + # + nodes = Node.([ Vec((-1.0, -1.0)), Vec((0.0, -1.0)), Vec((1.0, -1.0)), + Vec((-1.0, 0.0)), Vec((0.0, 0.0)), Vec((1.0, 0.0)), + Vec((-1.0, 1.0)), Vec((0.0, 1.0)), Vec((1.0, 1.0))]) + cells = [Triangle((1,2,4)), Triangle((2,5,4)), Quadrilateral((2,3,6,5)), Quadrilateral((4,5,8,7)), Quadrilateral((5,6,9,8))] + grid = Grid(cells, nodes) + addcellset!(grid, "bottom", OrderedSet((1,2,3))) + addcellset!(grid, "topleft", OrderedSet((4,))) + addcellset!(grid, "topright", OrderedSet((5,))) + + domain_names = ["bottom", "topleft", "topright"] + new_grid = insert_interfaces(grid, domain_names) + test_grid_data(grid, new_grid) +end + From 30167c843a31179d3e9a2eb7722f65e7fd6f3deb Mon Sep 17 00:00:00 2001 From: David Rollin Date: Tue, 7 Jul 2026 14:15:04 +0200 Subject: [PATCH 09/12] Improve tests and codecov --- src/cellvalues.jl | 18 ++++---- src/interpolations.jl | 1 + test/test_cellvalues.jl | 87 +++++++++++++++++++------------------ test/test_interpolations.jl | 5 +++ test/test_vtk_export.jl | 4 +- 5 files changed, 63 insertions(+), 52 deletions(-) diff --git a/src/cellvalues.jl b/src/cellvalues.jl index 9b37e2a..0b406bd 100755 --- a/src/cellvalues.jl +++ b/src/cellvalues.jl @@ -24,12 +24,12 @@ struct InterfaceCellValues{CV,TR,N} <: AbstractCellValues sides_and_baseindices::NTuple{N,Tuple{Symbol,Int}} R::TR # Union{AbstractVector,Nothing} Rotation matrix in quadrature points - function InterfaceCellValues(ip::IP, here::CV; use_same_cv::Bool, include_R::Val) where {IP<:InterfaceCellInterpolation, CV<:CellValues} + function InterfaceCellValues(ip::IP, here::CV; use_same_cv::Val, include_R::Val) where {IP<:InterfaceCellInterpolation, CV<:CellValues} N = getnbasefunctions(ip) sides_and_baseindices = Tuple( get_side_and_baseindex(ip, i) for i in 1:N ) base_indices_here = collect( get_interface_index(ip, :here, i) for i in 1:getnbasefunctions(ip.base) ) base_indices_there = collect( get_interface_index(ip, :there, i) for i in 1:getnbasefunctions(ip.base) ) - there = use_same_cv ? here : deepcopy(here) + there = use_same_cv === Val(true) ? here : deepcopy(here) R = if include_R === Val(false) nothing else @@ -38,13 +38,13 @@ struct InterfaceCellValues{CV,TR,N} <: AbstractCellValues end return new{CV,typeof(R),N}(here, there, base_indices_here, base_indices_there, sides_and_baseindices, R) end - function InterfaceCellValues(ip::IP, here::CV; use_same_cv, include_R::Val) where {IP<:VectorizedInterpolation{<:Any,<:Any,<:Any,<:InterfaceCellInterpolation}, CV<:CellValues} + function InterfaceCellValues(ip::IP, here::CV; use_same_cv::Val, include_R::Val) where {IP<:VectorizedInterpolation{<:Any,<:Any,<:Any,<:InterfaceCellInterpolation}, CV<:CellValues} N = getnbasefunctions(ip) sides_and_baseindices = Tuple( get_side_and_baseindex(ip, i) for i in 1:N ) ip = ip.ip base_indices_here = collect( get_interface_index(ip, :here, i) for i in 1:getnbasefunctions(ip.base) ) base_indices_there = collect( get_interface_index(ip, :there, i) for i in 1:getnbasefunctions(ip.base) ) - there = use_same_cv ? here : deepcopy(here) + there = use_same_cv === Val(true) ? here : deepcopy(here) R = if include_R === Val(false) nothing else @@ -60,17 +60,19 @@ InterfaceCellValues(qr::QuadratureRule, args...; kwargs...) = InterfaceCellValue function InterfaceCellValues(::Type{T}, qr::QuadratureRule, ip::InterfaceCellInterpolation, ip_geo::VectorizedInterpolation{sdim,<:Any,<:Any,<:InterfaceCellInterpolation} = default_geometric_interpolation(ip); - use_same_cv=true, include_R=Val(false), kwargs...) where {T, sdim} + use_same_cv=Val(true), include_R=Val(false), kwargs...) where {T, sdim} cv = CellValues(T, qr, ip.base, VectorizedInterpolation{sdim}(ip_geo.ip.base); kwargs...) - return InterfaceCellValues(ip, cv; use_same_cv=use_same_cv, include_R = isa(include_R, Bool) ? Val(include_R) : include_R) + return InterfaceCellValues(ip, cv; use_same_cv = isa(use_same_cv, Bool) ? Val(use_same_cv) : use_same_cv, + include_R = isa(include_R, Bool) ? Val(include_R) : include_R) end function InterfaceCellValues(::Type{T}, qr::QuadratureRule, ip::VectorizedInterpolation{vdim,<:Any,<:Any,<:InterfaceCellInterpolation}, ip_geo::VectorizedInterpolation{sdim,<:Any,<:Any,<:InterfaceCellInterpolation} = default_geometric_interpolation(ip); - use_same_cv=true, include_R=Val(false), kwargs...) where {T, vdim, sdim} + use_same_cv=Val(true), include_R=Val(false), kwargs...) where {T, vdim, sdim} cv = CellValues(T, qr, VectorizedInterpolation{vdim}(ip.ip.base), VectorizedInterpolation{sdim}(ip_geo.ip.base); kwargs...) - return InterfaceCellValues(ip, cv; use_same_cv=use_same_cv, include_R = isa(include_R, Bool) ? Val(include_R) : include_R) + return InterfaceCellValues(ip, cv; use_same_cv = isa(use_same_cv, Bool) ? Val(use_same_cv) : use_same_cv, + include_R = isa(include_R, Bool) ? Val(include_R) : include_R) end function InterfaceCellValues(::Type{T}, qr::QuadratureRule, diff --git a/src/interpolations.jl b/src/interpolations.jl index 07b7f9c..2701d01 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -93,6 +93,7 @@ function get_interface_index(ip::InterfaceCellInterpolation, side::Symbol, i::In nv = _nvertexdofs(ip.base) ne = _nedgedofs(ip.base) nf = _nfacedofs(ip.base) + @assert i > 0 if side == :here if i ≤ nv return i diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index fee05f0..3fd0245 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -5,47 +5,47 @@ (InterfaceCellInterpolation(Lagrange{RefTriangle, 1}()), InterfaceCellInterpolation(Lagrange{RefTriangle, 2}())), (InterfaceCellInterpolation(Lagrange{RefTriangle, 2}()), InterfaceCellInterpolation(Lagrange{RefTriangle, 2}())), ) - @test InterfaceCellValues(qr, fip, gip^3) isa InterfaceCellValues - @test InterfaceCellValues(qr, fip^3, gip^3) isa InterfaceCellValues - @test InterfaceCellValues(qr, fip^3, gip) isa InterfaceCellValues - @test InterfaceCellValues(qr, fip, gip) isa InterfaceCellValues - @test InterfaceCellValues(Float32, qr, fip) isa InterfaceCellValues - @test InterfaceCellValues(qr, fip) isa InterfaceCellValues - @test InterfaceCellValues(qr, fip; use_same_cv=false) isa InterfaceCellValues + @inferred InterfaceCellValues(qr, fip, gip^3) + #@inferred InterfaceCellValues(qr, fip^3, gip^3) + #@inferred InterfaceCellValues(qr, fip^3, gip) + @inferred InterfaceCellValues(qr, fip, gip) + @inferred InterfaceCellValues(Float32, qr, fip) + @inferred InterfaceCellValues(qr, fip) + @inferred InterfaceCellValues(qr, fip; use_same_cv=Val(false)) + @inferred InterfaceCellValues(qr, fip; use_same_cv=Val(true)) + @inferred InterfaceCellValues(qr, fip; include_R=Val(false)) + @inferred InterfaceCellValues(qr, fip; include_R=Val(true)) end ip = InterfaceCellInterpolation(Lagrange{RefTriangle, 1}()) - cv = InterfaceCellValues(qr, ip) + for cv in (InterfaceCellValues(qr, ip), InterfaceCellValues(qr, ip; use_same_cv=false), InterfaceCellValues(qr, ip; include_R=true)) + @test getnbasefunctions(cv) == 6 + @test Ferrite.getngeobasefunctions(cv) == 6 - @test getnbasefunctions(cv) == 6 - @test Ferrite.getngeobasefunctions(cv) == 6 - - x = repeat([rand(Vec{3}), rand(Vec{3}), rand(Vec{3})], 2) - reinit!(cv, x) - nbf = getnbasefunctions(cv) - here, there = rand(2) - u = vcat(ones(nbf÷2).*here, ones(nbf÷2).*there) - for qp in 1:getnquadpoints(cv) - @test function_value(cv, qp, u, true) ≈ here - @test function_value(cv, qp, u, false) ≈ there - @test all(abs.(function_gradient(cv, qp, u, true)) .≤ 1e-14) - @test all(abs.(function_gradient(cv, qp, u, false)) .≤ 1e-14) - @test function_value_average(cv, qp, u) ≈ (here + there)/2 - @test function_value_jump(cv, qp, u) ≈ there - here - @test all(abs.(function_gradient_average(cv, qp, u)) .≤ 1e-14) - @test all(abs.(function_gradient_jump(cv, qp, u)) .≤ 1e-14) - @test getdetJdV_average(cv, qp) == (getdetJdV(cv.here, qp) + getdetJdV(cv.there, qp)) / 2 - n = @allocated function_value_jump(cv, qp, u) - @test n == 0 + x = repeat([rand(Vec{3}), rand(Vec{3}), rand(Vec{3})], 2) + reinit!(cv, x) + nbf = getnbasefunctions(cv) + here, there = rand(2) + u = vcat(ones(nbf÷2).*here, ones(nbf÷2).*there) + for qp in 1:getnquadpoints(cv) + @test function_value(cv, qp, u, true) ≈ here + @test function_value(cv, qp, u, false) ≈ there + @test all(abs.(function_gradient(cv, qp, u, true)) .≤ 1e-14) + @test all(abs.(function_gradient(cv, qp, u, false)) .≤ 1e-14) + @test function_value_average(cv, qp, u) ≈ (here + there)/2 + @test function_value_jump(cv, qp, u) ≈ there - here + @test all(abs.(function_gradient_average(cv, qp, u)) .≤ 1e-14) + @test all(abs.(function_gradient_jump(cv, qp, u)) .≤ 1e-14) + @test getdetJdV_average(cv, qp) == (getdetJdV(cv.here, qp) + getdetJdV(cv.there, qp)) / 2 + n = @allocated function_value_jump(cv, qp, u) + @test n == 0 + end end qr = QuadratureRule{RefTriangle}(2) ip = InterfaceCellInterpolation(Lagrange{RefTriangle, 1}()) cv = InterfaceCellValues(qr, ip; include_R=true) - @inferred InterfaceCellValues(qr, ip) # default type stable - @inferred InterfaceCellValues(qr, ip; include_R=Val(false)) # no rotation type stable with Val - @inferred InterfaceCellValues(qr, ip; include_R=Val(true)) # with rotation type stable with Val. - + x = repeat([Vec{3}((1.0,0.0,0.0)), Vec{3}((0.0,1.0,0.0)), Vec{3}((0.0,0.0,1.0))], 2) n = Vec{3}(( sqrt(1/3), sqrt(1/3), sqrt(1/3))) t₁ = Vec{3}(( sqrt(1/2), 0.0, -sqrt(1/2))) @@ -60,16 +60,17 @@ end qr = QuadratureRule{RefLine}(2) - ip = InterfaceCellInterpolation(Lagrange{RefLine, 1}()) - cv = InterfaceCellValues(qr, ip; include_R=true) - x = repeat([Vec{2}((1.0,0.0)), Vec{2}((0.0,1.0))], 2) - n = Vec{2}((-sqrt(1/2), -sqrt(1/2))) - t = Vec{2}((-sqrt(1/2), sqrt(1/2))) - reinit!(cv, x) - for qp in 1:getnquadpoints(cv) - R = midplane_rotation(cv, qp) - @test tdot(R) ≈ one(R) # Fails e.g. when dx/dξ₁ not perpendicular to dx/ξ₂ - @test R⋅Vec{2}((1.0,0.0)) ≈ t - @test R⋅Vec{2}((0.0,1.0)) ≈ n + for ip in (InterfaceCellInterpolation(Lagrange{RefLine, 1}()), InterfaceCellInterpolation(Lagrange{RefLine, 1}())^2) + cv = InterfaceCellValues(qr, ip; include_R=true) + x = repeat([Vec{2}((1.0,0.0)), Vec{2}((0.0,1.0))], 2) + n = Vec{2}((-sqrt(1/2), -sqrt(1/2))) + t = Vec{2}((-sqrt(1/2), sqrt(1/2))) + reinit!(cv, x) + for qp in 1:getnquadpoints(cv) + R = midplane_rotation(cv, qp) + @test tdot(R) ≈ one(R) # Fails e.g. when dx/dξ₁ not perpendicular to dx/ξ₂ + @test R⋅Vec{2}((1.0,0.0)) ≈ t + @test R⋅Vec{2}((0.0,1.0)) ≈ n + end end end \ No newline at end of file diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 6e4e992..07648a3 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -26,4 +26,9 @@ expectedtype = InterfaceCellInterpolation{RefQuadrilateral, 1, Lagrange{RefLine,1}} @test Ferrite.default_geometric_interpolation(testcelltype) isa expectedtype @test Ferrite.default_geometric_interpolation(Ferrite.default_geometric_interpolation(testcelltype)) isa VectorizedInterpolation{2, RefQuadrilateral, <:Any, expectedtype} + + @test_throws AssertionError FerriteInterfaceElements.get_interface_index(ip, :here, 0) + @test_throws ArgumentError FerriteInterfaceElements.get_interface_index(ip, :here, 100) + @test_throws ArgumentError FerriteInterfaceElements.get_interface_index(ip, :there, 100) + @test_throws ArgumentError FerriteInterfaceElements.get_interface_index(ip, :test, 1) end diff --git a/test/test_vtk_export.jl b/test/test_vtk_export.jl index 3e9debd..b2f2f62 100644 --- a/test/test_vtk_export.jl +++ b/test/test_vtk_export.jl @@ -17,7 +17,9 @@ using Ferrite, FerriteInterfaceElements, OrderedCollections set_interface = getcellset(grid2, "interfaces") add!(SubDofHandler(dh, set_interface), :u, InterfaceCellInterpolation(Lagrange{FerriteInterfaceElements.getinterfaceshape(grid.cells[1]), 1}())) close!(dh) - VTKGridFile("debug.vtu", grid2) do vtk + temp = tempdir() + file = joinpath(temp, "output.vtu") + VTKGridFile(file, grid2) do vtk Ferrite.write_solution(vtk, dh, rand(ndofs(dh))) Ferrite.write_cellset(vtk, grid2) end From 65f85f484a213af655a21057222d009caa9b45a3 Mon Sep 17 00:00:00 2001 From: David Rollin Date: Tue, 7 Jul 2026 15:21:48 +0200 Subject: [PATCH 10/12] Support interfaces between pyramids and wedges --- Project.toml | 2 +- src/FerriteInterfaceElements.jl | 2 +- src/grid.jl | 32 +++++++++++++++++++++----------- test/test_cellvalues.jl | 1 + test/test_grid.jl | 19 +++++++++++++++++++ 5 files changed, 43 insertions(+), 13 deletions(-) diff --git a/Project.toml b/Project.toml index 5c45996..a52c2d9 100644 --- a/Project.toml +++ b/Project.toml @@ -11,8 +11,8 @@ VTKBase = "4004b06d-e244-455f-a6ce-a5f9919cc534" [compat] Ferrite = "1.4" OrderedCollections = "1" -julia = "1.12" VTKBase = "1.0.1" +julia = "1.12" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/FerriteInterfaceElements.jl b/src/FerriteInterfaceElements.jl index 7982bdc..78950c7 100644 --- a/src/FerriteInterfaceElements.jl +++ b/src/FerriteInterfaceElements.jl @@ -3,7 +3,7 @@ module FerriteInterfaceElements import Ferrite import Ferrite: AbstractCell, AbstractRefShape, AbstractCellValues, RefLine, RefQuadrilateral, RefTriangle, RefPrism, RefHexahedron, RefTetrahedron, - Line, QuadraticLine, Triangle, QuadraticTriangle, Quadrilateral, QuadraticQuadrilateral, Tetrahedron, Hexahedron, + Line, QuadraticLine, Triangle, QuadraticTriangle, Quadrilateral, QuadraticQuadrilateral, Tetrahedron, Hexahedron, Pyramid, Wedge, ScalarInterpolation, VectorizedInterpolation, Lagrange, CellValues, QuadratureRule, CellCache, Grid, ExclusiveTopology, FacetIndex, Vec, Tensor, MixedTensor2, getnbasefunctions, getngeobasefunctions, getorder, n_components, getrefshape, diff --git a/src/grid.jl b/src/grid.jl index ede67b8..cf140d2 100644 --- a/src/grid.jl +++ b/src/grid.jl @@ -97,7 +97,7 @@ end Return a suitable `InterfaceCell` connecting the facets with `nodes_here` and `nodes_there`. """ function create_interface_cell(::Type{C₁}, ::Type{C₂}, nodes_here, nodes_there) where {C₁,C₂} - Cbase = get_interface_base_cell_type(C₁, C₂) + Cbase = get_interface_base_cell_type(C₁, C₂, nodes_here) return InterfaceCell(Cbase(nodes_here), Cbase(nodes_there)) end @@ -107,14 +107,24 @@ end Return a suitable base type for connecting two cells of given type with an `InterfaceCell`. """ -get_interface_base_cell_type(::Type{Triangle}, ::Type{Triangle}) = Line -get_interface_base_cell_type(::Type{QuadraticTriangle}, ::Type{QuadraticTriangle}) = QuadraticLine -get_interface_base_cell_type(::Type{Quadrilateral}, ::Type{Quadrilateral}) = Line -get_interface_base_cell_type(::Type{QuadraticQuadrilateral}, ::Type{QuadraticQuadrilateral}) = QuadraticLine -get_interface_base_cell_type(::Type{Tetrahedron}, ::Type{Tetrahedron}) = Triangle -get_interface_base_cell_type(::Type{Hexahedron}, ::Type{Hexahedron}) = Quadrilateral +get_interface_base_cell_type(::Type{Triangle}, ::Type{Triangle}, ::Any) = Line +get_interface_base_cell_type(::Type{QuadraticTriangle}, ::Type{QuadraticTriangle}, ::Any) = QuadraticLine +get_interface_base_cell_type(::Type{Quadrilateral}, ::Type{Quadrilateral}, ::Any) = Line +get_interface_base_cell_type(::Type{QuadraticQuadrilateral}, ::Type{QuadraticQuadrilateral}, ::Any) = QuadraticLine +get_interface_base_cell_type(::Type{Tetrahedron}, ::Type{Tetrahedron}, ::Any) = Triangle +get_interface_base_cell_type(::Type{Hexahedron}, ::Type{Hexahedron}, ::Any) = Quadrilateral -get_interface_base_cell_type(::Type{Triangle}, ::Type{Quadrilateral}) = Line -get_interface_base_cell_type(::Type{Quadrilateral}, ::Type{Triangle}) = Line -get_interface_base_cell_type(::Type{QuadraticTriangle}, ::Type{QuadraticQuadrilateral}) = QuadraticLine -get_interface_base_cell_type(::Type{QuadraticQuadrilateral}, ::Type{QuadraticTriangle}) = QuadraticLine +get_interface_base_cell_type(::Type{Triangle}, ::Type{Quadrilateral}, ::Any) = Line +get_interface_base_cell_type(::Type{Quadrilateral}, ::Type{Triangle}, ::Any) = Line +get_interface_base_cell_type(::Type{QuadraticTriangle}, ::Type{QuadraticQuadrilateral}, ::Any) = QuadraticLine +get_interface_base_cell_type(::Type{QuadraticQuadrilateral}, ::Type{QuadraticTriangle}, ::Any) = QuadraticLine + +function get_interface_base_cell_type(::Type{C₁}, ::Type{C₂}, nodes) where {C₁<:Union{Pyramid,Wedge}, C₂<:Union{Pyramid,Wedge}} + if length(nodes) == 4 + return Quadrilateral + elseif length(nodes) == 3 + return Triangle + end + throw(ErrorException("No feasible base type for InterfaceCell!")) + return nothing +end diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 3fd0245..c23d272 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -37,6 +37,7 @@ @test all(abs.(function_gradient_average(cv, qp, u)) .≤ 1e-14) @test all(abs.(function_gradient_jump(cv, qp, u)) .≤ 1e-14) @test getdetJdV_average(cv, qp) == (getdetJdV(cv.here, qp) + getdetJdV(cv.there, qp)) / 2 + precompile(function_value_jump, (typeof(cv), typeof(qp), typeof(u))) n = @allocated function_value_jump(cv, qp, u) @test n == 0 end diff --git a/test/test_grid.jl b/test/test_grid.jl index 6015711..b1fc2b6 100644 --- a/test/test_grid.jl +++ b/test/test_grid.jl @@ -161,3 +161,22 @@ end test_grid_data(grid, new_grid) end +@testset "Inserting interfaces in 3D for a mixed grid" begin + nodes = Node.([ Vec((-1.0, -1.0, -1.0)), Vec(( 0.0, -1.0, -1.0)), Vec(( 0.0, 0.0, -1.0)), Vec((-1.0, 0.0, -1.0)), # bottom of pyramids 1+2 + Vec((-1.0, -1.0, 0.0)), # tip of pyramid 1 + Vec((-1.0, -1.0, -2.0)), # tip of pyramid 2 + Vec((-1.0, -2.0, -1.0)), Vec(( 0.0, -2.0, -1.0)), # bottom of pyramid 3 + Vec((-1.0, -1.0, -2.0)), Vec((-1.0, -2.0, -2.0)), # wedge 1 + ]) + cells = [Pyramid((1,2,3,4,5)), Pyramid((1,2,3,4,6)), Pyramid((1,2,8,7,5)), Wedge((1,2,9,8,7,10))] + grid = Grid(cells, nodes) + addcellset!(grid, "p 1", OrderedSet((1,))) + addcellset!(grid, "p 2", OrderedSet((2,))) + addcellset!(grid, "p 3", OrderedSet((3,))) + addcellset!(grid, "w 1", OrderedSet((4,))) + + domain_names = ["p 1", "p 2", "p 3", "w 1"] + new_grid = insert_interfaces(grid, domain_names) + @test length(new_grid.cells) == 7 + @test length(new_grid.nodes) == 21 +end \ No newline at end of file From 428c68df600ad38c83aec92619efc5800cb73036 Mon Sep 17 00:00:00 2001 From: David Rollin Date: Wed, 8 Jul 2026 11:10:44 +0200 Subject: [PATCH 11/12] Extended interface insertion to defining specific interfaces --- src/grid.jl | 158 ++++++++++++++++++++++++++++------------------ test/test_grid.jl | 27 ++++++-- 2 files changed, 117 insertions(+), 68 deletions(-) diff --git a/src/grid.jl b/src/grid.jl index cf140d2..427b586 100644 --- a/src/grid.jl +++ b/src/grid.jl @@ -2,70 +2,81 @@ # Inserting cells into a grid ###################################################################### """ - insert_interfaces(grid, domain_names; topology=ExclusiveTopology(grid)) + insert_interfaces(grid::Grid, domain_names::Vector{String}; kwargs...) + insert_interfaces(grid::Grid, interfaces::Dict{String,Tuple{String,String}}; kwargs...) + +Return a new grid with `InterfaceCell`s inserted according to `domain_names` or `interfaces`. +The new grid provides additional cell sets. The set `"interfaces"` contains all new `InterfaceCell`s. -Return a new grid with `InterfaceCell`s inserted between the domains defined by `domain_names`. -The new grid provides additional cell sets. The set `"interfaces"` contains all new `InterfaceCell`s -and two sets are provided for each combination of domain names: `"domain1-domain2-interface"` and -`"domain2-domain1-interface"` both using the same `Set`. +When using `domain_names`, for each combination of the corresponding cellsests interfaces will be inserted +and a new cellset will be provided, which can be accessed using both of the following names: +`"domain1-domain2-interface"` and `"domain2-domain1-interface"`. + +When using `interfaces`, each value of the `Dict` defines a pair of cellsets between which interfaces +will be inserted and collected in a new cellset which can be accessed by the corresponding key. """ -function insert_interfaces(grid, domain_names; topology=ExclusiveTopology(grid)) - cellsets = Dict(name => getcellset(grid, name) for name in domain_names) - pairs = OrderedSet{Pair{String,OrderedSet{Int}}}() # prepare a cellset for each interface between domains - for i in eachindex(domain_names) - for j in i+1:length(domain_names) - set = OrderedSet{Int}() - # add the set with two names to allow for both orders of domain names - push!(pairs, "$(domain_names[i])-$(domain_names[j])-interface" => set) - push!(pairs, "$(domain_names[j])-$(domain_names[i])-interface" => set) - end - end - interfacesets = Dict(pairs) - node_mapping = Dict(name => Dict{Int, Int}() for name in domain_names) +function insert_interfaces(grid::Grid, interfaces::Dict{String,Tuple{String,String}}; kwargs...) + interfaces = [(name, domains[1], domains[2]) for (name, domains) in pairs(interfaces)] + return _insert_interfaces(grid, interfaces; kwargs...) +end +function insert_interfaces(grid::Grid, domain_names::Vector{String}; kwargs...) + interfaces = [begin + names = ("$(domain_names[i])-$(domain_names[j])-interface", "$(domain_names[j])-$(domain_names[i])-interface") + (names, domain_names[i], domain_names[j]) + end for i in 1:length(domain_names) for j in i+1:length(domain_names)] + return _insert_interfaces(grid, interfaces; kwargs...) +end + + + +function _insert_interfaces(grid::Grid, interfaces::Vector{Tuple{T,String,String}}; topology=ExclusiveTopology(grid)) where {T<:Union{String,Tuple{String,String}}} + relevant_domains, cellsets = _prepare_cellsets(grid, interfaces) + interfacesets = _prepare_interfacesets(interfaces) + node_mapping = Dict(name => Dict{Int, Int}() for name in keys(cellsets)) nodes = copy(grid.nodes) cells_generic = Vector{Ferrite.AbstractCell}(grid.cells) # copies - while ! isempty(cellsets) - name, cellset = pop!(cellsets) - for cellid in cellset - cell = getcells(grid, cellid) - for facetid in 1:length(facets(cell)) - facet_neighbors = getneighborhood(topology, grid, FacetIndex(cellid, facetid)) + for (name, domain_h, domain_t) in interfaces + cellset_h = cellsets[domain_h] + cellset_t = cellsets[domain_t] + for cellid_h in cellset_h + cell_h = getcells(grid, cellid_h) + for facetid_h in 1:length(facets(cell_h)) + facet_neighbors = getneighborhood(topology, grid, FacetIndex(cellid_h, facetid_h)) isempty(facet_neighbors) && continue - (_cellid, _facetid) = only(facet_neighbors) # should only ever be one neighboring face - for (_name, _cellset) in cellsets # only "other" cellsets left - if _cellid ∈ _cellset # interface detected - facetnodeids = Ferrite.facets(cell)[facetid] # original nodeids - for nodeid in facetnodeids - new_nodeid = get(node_mapping[name], nodeid, nothing) - _new_nodeid = get(node_mapping[_name], nodeid, nothing) - # generate missing duplicate nodes - if isnothing(new_nodeid) && isnothing(_new_nodeid) - # decide for main side, genererate new node - node_mapping[name][nodeid] = nodeid # main node - # new node - push!(nodes, nodes[nodeid]) - node_mapping[_name][nodeid] = length(nodes) # main node - elseif !isnothing(new_nodeid) && isnothing(_new_nodeid) - # node has been duplicated at least once before, so it already has a main - # generate a new node for (_name, _cellset) - push!(nodes, nodes[nodeid]) - node_mapping[_name][nodeid] = length(nodes) # main node - elseif isnothing(new_nodeid) && ! isnothing(_new_nodeid) - # node has been duplicated at least once, so it already has a main grain - # generate a new node for (name, cellset) - push!(nodes, nodes[nodeid]) - node_mapping[name][nodeid] = length(nodes) # main node - end + (cellid_t, facetid_t) = only(facet_neighbors) # should only ever be one neighboring face + if cellid_t in cellset_t # relevant interface detected + facetnodeids = Ferrite.facets(cell_h)[facetid_h] # original nodeids + for nodeid in facetnodeids + new_nodeid_h = get(node_mapping[domain_h], nodeid, nothing) + new_nodeid_t = get(node_mapping[domain_t], nodeid, nothing) + # generate missing duplicate nodes + if isnothing(new_nodeid_h) && isnothing(new_nodeid_t) + # decide for main side, genererate new node + node_mapping[domain_h][nodeid] = nodeid # main node + # new node + push!(nodes, nodes[nodeid]) + node_mapping[domain_t][nodeid] = length(nodes) # main node + elseif !isnothing(new_nodeid_h) && isnothing(new_nodeid_t) + # node has been duplicated at least once before, so it already has a main + # generate a new node for (_name, _cellset) + push!(nodes, nodes[nodeid]) + node_mapping[domain_t][nodeid] = length(nodes) # main node + elseif isnothing(new_nodeid_h) && ! isnothing(new_nodeid_t) + # node has been duplicated at least once, so it already has a main grain + # generate a new node for (name, cellset) + push!(nodes, nodes[nodeid]) + node_mapping[domain_h][nodeid] = length(nodes) # main node end - new_nodeids = Tuple(node_mapping[name][i] for i in facetnodeids) - _new_nodeids = Tuple(node_mapping[_name][i] for i in facetnodeids) - # generate new cell - interface_cell = create_interface_cell(typeof(cell), typeof(getcells(grid, _cellid)), new_nodeids, _new_nodeids) - push!(cells_generic, interface_cell) - push!(interfacesets["$(name)-$(_name)-interface"], length(cells_generic)) end + new_nodeids_h = Tuple(node_mapping[domain_h][i] for i in facetnodeids) + new_nodeids_t = Tuple(node_mapping[domain_t][i] for i in facetnodeids) + # generate new cell + cell_t = getcells(grid, cellid_t) + interface_cell = create_interface_cell(typeof(cell_h), typeof(cell_t), new_nodeids_h, new_nodeids_t) + push!(cells_generic, interface_cell) + _add_interfacecell!(interfacesets, length(cells_generic), name) end end end @@ -76,11 +87,11 @@ function insert_interfaces(grid, domain_names; topology=ExclusiveTopology(grid)) cells = convert(Array{cell_type}, cells_generic) # adjust original cells to new node numbering - for name in domain_names - cellset = getcellset(grid, name) + for domain in relevant_domains + cellset = getcellset(grid, domain) for cellid in cellset cell = getcells(grid, cellid) - cells[cellid] = typeof(cell)(map(n->get(node_mapping[name], n, n), cell.nodes)) + cells[cellid] = typeof(cell)(map(n -> get(node_mapping[domain], n, n), cell.nodes)) end end @@ -91,14 +102,35 @@ function insert_interfaces(grid, domain_names; topology=ExclusiveTopology(grid)) return new_grid end +function _prepare_cellsets(grid::Grid, interfaces::Vector{Tuple{T,String,String}}) where {T<:Union{String,Tuple{String,String}}} + relevant_domains = Set{String}() + for (_, domain_h, domain_t) in interfaces + push!(relevant_domains, domain_h) + push!(relevant_domains, domain_t) + end + return relevant_domains, Dict(name => getcellset(grid, name) for name in relevant_domains) +end + +function _prepare_interfacesets(interfaces::Vector{Tuple{String,String,String}}) + return Dict([ name => OrderedSet{Int}() for (name,_,_) in interfaces ]) +end +function _prepare_interfacesets(interfaces::Vector{Tuple{Tuple{String,String},String,String}}) + sets = [OrderedSet{Int}() for _ in interfaces] + return Dict([ name => set for (set, (names,_,_)) in zip(sets, interfaces) for name in names ]) +end + +_add_interfacecell!(interfacesets::Dict{String,OrderedSet{Int}}, cellid::Int, name::String) = push!(interfacesets[name], cellid) +_add_interfacecell!(interfacesets::Dict{String,OrderedSet{Int}}, cellid::Int, name::Tuple{String,String}) = push!(interfacesets[name[1]], cellid) + + """ - create_interface_cell(::Type{C₁}, ::Type{C₂}, nodes_here, nodes_there) where {C₁,C₂} + create_interface_cell(::Type{C₁}, ::Type{C₂}, nodes_h, nodes_t) where {C₁,C₂} -Return a suitable `InterfaceCell` connecting the facets with `nodes_here` and `nodes_there`. +Return a suitable `InterfaceCell` connecting the facets with `nodes_h` and `nodes_t`. """ -function create_interface_cell(::Type{C₁}, ::Type{C₂}, nodes_here, nodes_there) where {C₁,C₂} - Cbase = get_interface_base_cell_type(C₁, C₂, nodes_here) - return InterfaceCell(Cbase(nodes_here), Cbase(nodes_there)) +function create_interface_cell(::Type{C₁}, ::Type{C₂}, nodes_h, nodes_t) where {C₁,C₂} + Cbase = get_interface_base_cell_type(C₁, C₂, nodes_h) + return InterfaceCell(Cbase(nodes_h), Cbase(nodes_t)) end diff --git a/test/test_grid.jl b/test/test_grid.jl index b1fc2b6..4cc060b 100644 --- a/test/test_grid.jl +++ b/test/test_grid.jl @@ -85,7 +85,7 @@ end # |\ c6 |\ c8 | |\ c6 | c|\ c8 | # | \ | \ | | \ |11| \ | # | c5 \| c7 \| | c5 \| | c7 \| - # 4 ___ 5 ___ 6 ---> 11___10__13___12 + # 4 ___ 5 ___ 6 ---> 11___10__13___12 # |\ c2 |\ c4 | | c9 \/ c10 | # | \ | \ | 4 ___ 5 ___ 6 # | c1 \| c3 \| |\ c2 | \ c4 | @@ -101,6 +101,23 @@ end domain_names = ["bottom", "topleft", "topright"] new_grid = insert_interfaces(grid, domain_names) test_grid_data(grid, new_grid) + + # 7 ___ 8 ___ 9 7 ___ 8__13___ 9 + # |\ c6 |\ c8 | |\ c6 | c|\ c8 | + # | \ | \ | | \ |10| \ | + # | c5 \| c7 \| | c5 \| | c7 \| + # 4 ___ 5 ___ 6 ---> 11___10__12___ | + # |\ c2 |\ c4 | | c9 \/ 6 + # | \ | \ | 4 ___ 5 ____/| + # | c1 \| c3 \| |\ c2 | \ c4 | + # 1 ___ 2 ___ 3 | \ | \ | + # | c1 \ | c3 \| + # 1 ____ 2 ____ 3 + # + new_grid = insert_interfaces(grid, Dict(["A" => ("bottom", "topleft"), "B" => ("topright", "topleft")])) + @test length(new_grid.nodes) == 13 + @test length(new_grid.cells) == 10 + @test length(new_grid.cellsets) == 6 end @testset "Inserting interfaces in 3D" begin @@ -128,10 +145,10 @@ end setdiff!(set1, set2, set3) setdiff!(set3, set2) - newgrid = insert_interfaces(grid, ["1", "2", "3"]) - @test length(getcells(newgrid, "1-2-interface")) == 384 - @test length(getcells(newgrid, "1-3-interface")) == 192 - @test length(getcells(newgrid, "2-3-interface")) == 48 + new_grid = insert_interfaces(grid, ["1", "2", "3"]) + @test length(getcells(new_grid, "1-2-interface")) == 384 + @test length(getcells(new_grid, "1-3-interface")) == 192 + @test length(getcells(new_grid, "2-3-interface")) == 48 end @testset "Inserting interfaces in 2D for a mixed grid" begin From 4a79dfc7d26bd03f108a8efec6c1465d69c6b985 Mon Sep 17 00:00:00 2001 From: David Rollin Date: Thu, 9 Jul 2026 09:10:49 +0200 Subject: [PATCH 12/12] Fix test for julia 1.10 --- test/test_cellvalues.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index c23d272..ee10a5c 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -18,7 +18,7 @@ end ip = InterfaceCellInterpolation(Lagrange{RefTriangle, 1}()) - for cv in (InterfaceCellValues(qr, ip), InterfaceCellValues(qr, ip; use_same_cv=false), InterfaceCellValues(qr, ip; include_R=true)) + for cv in (InterfaceCellValues(qr, ip), InterfaceCellValues(qr, ip; use_same_cv=Val(false)), InterfaceCellValues(qr, ip; include_R=Val(true))) @test getnbasefunctions(cv) == 6 @test Ferrite.getngeobasefunctions(cv) == 6 @@ -37,7 +37,6 @@ @test all(abs.(function_gradient_average(cv, qp, u)) .≤ 1e-14) @test all(abs.(function_gradient_jump(cv, qp, u)) .≤ 1e-14) @test getdetJdV_average(cv, qp) == (getdetJdV(cv.here, qp) + getdetJdV(cv.there, qp)) / 2 - precompile(function_value_jump, (typeof(cv), typeof(qp), typeof(u))) n = @allocated function_value_jump(cv, qp, u) @test n == 0 end