From d21111ceb57a40386421e57e2ac80a4d89a34924 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Sat, 14 Oct 2023 14:14:09 +0200 Subject: [PATCH 1/7] Add docstring to SubDofHandler and improve docstring for DofHandler (#825) --- docs/src/reference/dofhandler.md | 1 + src/Dofs/DofHandler.jl | 60 +++++++++++++++++++++++++++----- 2 files changed, 52 insertions(+), 9 deletions(-) diff --git a/docs/src/reference/dofhandler.md b/docs/src/reference/dofhandler.md index 0bd5a839f5..f93deb2810 100644 --- a/docs/src/reference/dofhandler.md +++ b/docs/src/reference/dofhandler.md @@ -12,6 +12,7 @@ SubDofHandler ## Adding fields to the DofHandlers ```@docs add!(::DofHandler, ::Symbol, ::Interpolation) +add!(::SubDofHandler, ::Symbol, ::Interpolation) close!(::DofHandler) ``` diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index a653c8a563..8a0087f2d7 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -12,7 +12,6 @@ Access some grid representation for the dof handler. """ get_grid(dh::AbstractDofHandler) - struct SubDofHandler{DH} <: AbstractDofHandler # From constructor dh::DH @@ -26,8 +25,36 @@ struct SubDofHandler{DH} <: AbstractDofHandler # const dof_ranges::Vector{UnitRange{Int}} # TODO: Why not? end -# TODO: Should be an inner constructor. +""" + SubDofHandler(dh::AbstractDofHandler, cellset::Set{Int}) + +Create an `sdh::SubDofHandler` from the parent `dh`, pertaining to the +cells in `cellset`. This allows you to add fields to parts of the domain, or using +different interpolations or cell types (e.g. `Triangles` and `Quadrilaterals`). All +fields and cell types must be the same in one `SubDofHandler`. + +After construction any number of discrete fields can be added to the SubDofHandler using +[`add!`](@ref). Construction is finalized by calling [`close!`](@ref) on the parent `dh`. + +# Examples +We assume we have a `grid` containing "Triangle" and "Quadrilateral" cells, +including the cellsets "triangles" and "quadilaterals" for to these cells. +```julia +dh = DofHandler(grid) + +sdh_tri = SubDofHandler(dh, getcellset(grid, "triangles")) +ip_tri = Lagrange{RefTriangle, 2}()^2 # vector interpolation for a field u +add!(sdh_tri, :u, ip_tri) + +sdh_quad = SubDofHandler(dh, getcellset(grid, "quadilaterals")) +ip_quad = Lagrange{RefQuadrilateral, 2}()^2 # vector interpolation for a field u +add!(sdh_quad, :u, ip_quad) + +close!(dh) # Finalize by closing the parent +``` +""" function SubDofHandler(dh::DH, cellset) where {DH <: AbstractDofHandler} + # TODO: Should be an inner constructor. isclosed(dh) && error("DofHandler already closed") # Compute the celltype and make sure all elements have the same one CT = getcelltype(dh.grid, first(cellset)) @@ -66,13 +93,6 @@ function _print_field_information(io::IO, mime::MIME"text/plain", sdh::SubDofHan end end -""" - DofHandler(grid::Grid) - -Construct a `DofHandler` based on `grid`. Supports: -- `Grid`s with or without concrete element type (E.g. "mixed" grids with several different element types.) -- One or several fields, which can live on the whole domain or on subsets of the `Grid`. -""" struct DofHandler{dim,G<:AbstractGrid{dim}} <: AbstractDofHandler subdofhandlers::Vector{SubDofHandler{DofHandler{dim, G}}} field_names::Vector{Symbol} @@ -86,6 +106,28 @@ struct DofHandler{dim,G<:AbstractGrid{dim}} <: AbstractDofHandler ndofs::ScalarWrapper{Int} end +""" + DofHandler(grid::Grid) + +Construct a `DofHandler` based on the grid `grid`. + +After construction any number of discrete fields can be added to the DofHandler using +[`add!`](@ref). Construction is finalized by calling [`close!`](@ref). + +By default fields are added to all elements of the grid. Refer to [`SubDofHandler`](@ref) +for restricting fields to subdomains of the grid. + +# Examples + +```julia +dh = DofHandler(grid) +ip_u = Lagrange{RefTriangle, 2}()^2 # vector interpolation for a field u +ip_p = Lagrange{RefTriangle, 1}() # scalar interpolation for a field p +add!(dh, :u, ip_u) +add!(dh, :p, ip_p) +close!(dh) +``` +""" function DofHandler(grid::G) where {dim, G <: AbstractGrid{dim}} ncells = getncells(grid) sdhs = SubDofHandler{DofHandler{dim, G}}[] From fe44e7f6297158e9bbb62fd4bc9adaa1d79929f2 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Mon, 16 Oct 2023 09:43:51 +0200 Subject: [PATCH 2/7] Batch evaluation of shape function values/gradients (#819) --- docs/src/devdocs/interpolations.md | 3 +++ src/FEValues/cell_values.jl | 8 ++---- src/FEValues/face_values.jl | 14 ++++------- src/PointEval/PointEvalHandler.jl | 1 + src/PointEval/point_values.jl | 9 ++----- src/interpolations.jl | 40 ++++++++++++++++++++++++++++++ test/test_cellvalues.jl | 2 +- 7 files changed, 54 insertions(+), 23 deletions(-) diff --git a/docs/src/devdocs/interpolations.md b/docs/src/devdocs/interpolations.md index 23b3ec1a0a..b1016b054c 100644 --- a/docs/src/devdocs/interpolations.md +++ b/docs/src/devdocs/interpolations.md @@ -15,6 +15,9 @@ Ferrite.shape_gradient(::Interpolation, ::Vec, ::Int) Ferrite.shape_gradient_and_value Ferrite.boundarydof_indices Ferrite.dirichlet_boundarydof_indices +Ferrite.shape_values! +Ferrite.shape_gradients! +Ferrite.shape_gradients_and_values! ``` ### Required methods to implement for all subtypes of `Interpolation` to define a new finite element diff --git a/src/FEValues/cell_values.jl b/src/FEValues/cell_values.jl index 5ef93b8cff..e5d33f0aee 100644 --- a/src/FEValues/cell_values.jl +++ b/src/FEValues/cell_values.jl @@ -68,12 +68,8 @@ function CellValues{IP, N_t, dNdx_t, dNdξ_t, T, dMdξ_t, QR, GIP}(qr::QR, ip::I dMdξ = fill(zero(dMdξ_t) * T(NaN), n_geom_basefuncs, n_qpoints) for (qp, ξ) in pairs(getpoints(qr)) - for basefunc in 1:n_func_basefuncs - dNdξ[basefunc, qp], N[basefunc, qp] = shape_gradient_and_value(ip, ξ, basefunc) - end - for basefunc in 1:n_geom_basefuncs - dMdξ[basefunc, qp], M[basefunc, qp] = shape_gradient_and_value(gip, ξ, basefunc) - end + shape_gradients_and_values!(@view(dNdξ[:, qp]), @view(N[:, qp]), ip, ξ) + shape_gradients_and_values!(@view(dMdξ[:, qp]), @view(M[:, qp]), gip, ξ) end detJdV = fill(T(NaN), n_qpoints) diff --git a/src/FEValues/face_values.jl b/src/FEValues/face_values.jl index 5ba510330f..e520d8bce1 100644 --- a/src/FEValues/face_values.jl +++ b/src/FEValues/face_values.jl @@ -71,12 +71,8 @@ function FaceValues{IP, N_t, dNdx_t, dNdξ_t, T, dMdξ_t, QR, Normal_t, GIP}(qr: dMdξ = fill(zero(dMdξ_t) * T(NaN), n_geom_basefuncs, max_n_qpoints, n_faces) for face in 1:n_faces, (qp, ξ) in pairs(getpoints(qr, face)) - for basefunc in 1:n_func_basefuncs - dNdξ[basefunc, qp, face], N[basefunc, qp, face] = shape_gradient_and_value(ip, ξ, basefunc) - end - for basefunc in 1:n_geom_basefuncs - dMdξ[basefunc, qp, face], M[basefunc, qp, face] = shape_gradient_and_value(gip, ξ, basefunc) - end + shape_gradients_and_values!(@view(dNdξ[:, qp, face]), @view(N[:, qp, face]), ip, ξ) + shape_gradients_and_values!(@view(dMdξ[:, qp, face]), @view(M[:, qp, face]), gip, ξ) end detJdV = fill(T(NaN), max_n_qpoints, n_faces) @@ -223,8 +219,8 @@ function BCValues(::Type{T}, func_interpol::Interpolation{refshape}, geom_interp nqp = zeros(Int,n_boundary_entities) for n_boundary_entity in 1:n_boundary_entities - for (qp, ξ) in enumerate(qrs[n_boundary_entity].points), i in 1:n_geom_basefuncs - M[i, qp, n_boundary_entity] = shape_value(geom_interpol, ξ, i) + for (qp, ξ) in pairs(qrs[n_boundary_entity].points) + shape_values!(@view(M[:, qp, n_boundary_entity]), geom_interpol, ξ) end nqp[n_boundary_entity] = length(qrs[n_boundary_entity].points) end @@ -262,4 +258,4 @@ function Base.show(io::IO, m::MIME"text/plain", fv::FaceValues) print(io, "- Function interpolation: "); show(io, m, fv.func_interp) println(io) print(io, "- Geometric interpolation: "); show(io, m, fv.geo_interp) -end \ No newline at end of file +end diff --git a/src/PointEval/PointEvalHandler.jl b/src/PointEval/PointEvalHandler.jl index 62b2dfdb49..db9d16fc93 100644 --- a/src/PointEval/PointEvalHandler.jl +++ b/src/PointEval/PointEvalHandler.jl @@ -125,6 +125,7 @@ function find_local_coordinate(interpolation, cell_coordinates::Vector{V}, globa for _ in 1:max_iters global_guess = zero(V) J = zero(Tensor{2, dim, T}) + # TODO batched eval after 764 is merged. for j in 1:n_basefuncs dNdξ, N = shape_gradient_and_value(interpolation, local_guess, j) global_guess += N * cell_coordinates[j] diff --git a/src/PointEval/point_values.jl b/src/PointEval/point_values.jl index a49b57ce76..15873a5bbb 100644 --- a/src/PointEval/point_values.jl +++ b/src/PointEval/point_values.jl @@ -59,9 +59,7 @@ function_symmetric_gradient(pv::PointValues, u::AbstractVector, args...) = function reinit!(pv::PointValues, x::AbstractVector{<:Vec{D}}, ξ::Vec{D}) where {D} qp = 1 # PointValues only have a single qp # TODO: Does M need to be updated too? - for i in 1:getnbasefunctions(pv.cv.ip) - pv.cv.dNdξ[i, qp], pv.cv.N[i, qp] = shape_gradient_and_value(pv.cv.ip, ξ, i) - end + shape_gradients_and_values!(@view(pv.cv.dNdξ[:, qp]), @view(pv.cv.N[:, qp]), pv.cv.ip, ξ) reinit!(pv.cv, x) return nothing end @@ -86,9 +84,6 @@ shape_value(pv::PointValuesInternal, qp::Int, i::Int) = (@assert qp == 1; pv.N[i # allow on-the-fly updating function reinit!(pv::PointValuesInternal{IP}, coord::Vec{dim}) where {dim, shape <: AbstractRefShape{dim}, IP <: Interpolation{shape}} - n_func_basefuncs = getnbasefunctions(pv.ip) - for i in 1:n_func_basefuncs - pv.N[i] = shape_value(pv.ip, coord, i) - end + shape_values!(pv.N, pv.ip, coord) return nothing end diff --git a/src/interpolations.jl b/src/interpolations.jl index 649107c177..4e1e70bf11 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -204,6 +204,46 @@ getnbasefunctions(::Interpolation) # edgedof: dof on a line between 2 vertices (i.e. "corners") (3D only) # celldof: dof that is local to the element +""" + shape_values!(values::AbstractArray{T}, ip::Interpolation, ξ::Vec) + +Evaluate all shape functions of `ip` at once at the reference point `ξ` and store them in +`values`. +""" +@propagate_inbounds function shape_values!(values::AT, ip::IP, ξ::Vec) where {IP <: Interpolation, AT <: AbstractArray} + @boundscheck checkbounds(values, 1:getnbasefunctions(ip)) + @inbounds for i in 1:getnbasefunctions(ip) + values[i] = shape_value(ip, ξ, i) + end +end + +""" + shape_gradients!(gradients::AbstractArray, ip::Interpolation, ξ::Vec) + +Evaluate all shape function gradients of `ip` at once at the reference point `ξ` and store +them in `gradients`. +""" +function shape_gradients!(gradients::AT, ip::IP, ξ::Vec) where {IP <: Interpolation, AT <: AbstractArray} + @boundscheck checkbounds(gradients, 1:getnbasefunctions(ip)) + @inbounds for i in 1:getnbasefunctions(ip) + gradients[i] = shape_gradient(ip, ξ, i) + end +end + +""" + shape_gradients_and_values!(gradients::AbstractArray, values::AbstractArray, ip::Interpolation, ξ::Vec) + +Evaluate all shape function gradients and values of `ip` at once at the reference point `ξ` +and store them in `values`. +""" +function shape_gradients_and_values!(gradients::GAT, values::SAT, ip::IP, ξ::Vec) where {IP <: Interpolation, SAT <: AbstractArray, GAT <: AbstractArray} + @boundscheck checkbounds(gradients, 1:getnbasefunctions(ip)) + @boundscheck checkbounds(values, 1:getnbasefunctions(ip)) + @inbounds for i in 1:getnbasefunctions(ip) + gradients[i], values[i] = shape_gradient_and_value(ip, ξ, i) + end +end + """ shape_value(ip::Interpolation, ξ::Vec, i::Int) diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 714741ed05..7db59cc53b 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -1,5 +1,5 @@ @testset "CellValues" begin -for (scalar_interpol, quad_rule) in ( +@testset "ip=$scalar_interpol quad_rule=$(typeof(quad_rule))" for (scalar_interpol, quad_rule) in ( (Lagrange{RefLine, 1}(), QuadratureRule{RefLine}(2)), (Lagrange{RefLine, 2}(), QuadratureRule{RefLine}(2)), (Lagrange{RefQuadrilateral, 1}(), QuadratureRule{RefQuadrilateral}(2)), From 9e05989647dd4d18589490447676c1d8d14f8092 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Mon, 6 Nov 2023 22:32:26 +0100 Subject: [PATCH 3/7] Uniform point sampling for reference geometry and interpolation tests, fixes #811 (#818) --- test/test_interpolations.jl | 54 +++++++++++++++++++++---------------- test/test_utils.jl | 22 +++++++++++++++ 2 files changed, 53 insertions(+), 23 deletions(-) diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index b6653b4491..ee488cbc69 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -37,7 +37,6 @@ # CrouzeixRaviart{RefTriangle, 1}(), ) - # Test of utility functions ref_dim = Ferrite.getdim(interpolation) ref_shape = Ferrite.getrefshape(interpolation) @@ -47,19 +46,28 @@ # Note that not every element formulation exists for every order and dimension. applicable(Ferrite.getlowerorder, interpolation) && @test typeof(Ferrite.getlowerorder(interpolation)) <: Interpolation{ref_shape,func_order-1} - # Check partition of unity at random point. - n_basefuncs = getnbasefunctions(interpolation) - x = rand(Tensor{1, ref_dim}) - f = (x) -> [shape_value(interpolation, Tensor{1, ref_dim}(x), i) for i in 1:n_basefuncs] - @test vec(ForwardDiff.jacobian(f, Array(x))') ≈ - reinterpret(Float64, [shape_gradient(interpolation, x, i) for i in 1:n_basefuncs]) - @test sum([shape_value(interpolation, x, i) for i in 1:n_basefuncs]) ≈ 1.0 + n_basefuncs = getnbasefunctions(interpolation) + coords = Ferrite.reference_coordinates(interpolation) + @test length(coords) == n_basefuncs + f(x) = [shape_value(interpolation, Tensor{1, ref_dim}(x), i) for i in 1:n_basefuncs] - # Check if the important functions are consistent - coords = Ferrite.reference_coordinates(interpolation) - @test length(coords) == n_basefuncs - @test shape_value(interpolation, x, n_basefuncs) == shape_value(interpolation, x, n_basefuncs) - @test_throws ArgumentError shape_value(interpolation, x, n_basefuncs+1) + #TODO prefer this test style after 1.6 is removed from CI + # @testset let x = sample_random_point(ref_shape) # not compatible with Julia 1.6 + x = sample_random_point(ref_shape) + random_point_testset = @testset "Random point test" begin + # Check gradient evaluation + @test vec(ForwardDiff.jacobian(f, Array(x))') ≈ + reinterpret(Float64, [shape_gradient(interpolation, x, i) for i in 1:n_basefuncs]) + # Check partition of unity at random point. + @test sum([shape_value(interpolation, x, i) for i in 1:n_basefuncs]) ≈ 1.0 + # Check if the important functions are consistent + @test_throws ArgumentError shape_value(interpolation, x, n_basefuncs+1) + # Idempotency test + @test shape_value(interpolation, x, n_basefuncs) == shape_value(interpolation, x, n_basefuncs) + end + # Remove after 1.6 is removed from CI (see above) + # Show coordinate in case failure (see issue #811) + !isempty(random_point_testset.results) && println("^^^^^Random point test failed at $x for $interpolation !^^^^^") # Test whether we have for each entity corresponding dof indices (possibly empty) @test length(Ferrite.vertexdof_indices(interpolation)) == Ferrite.nvertices(interpolation) @@ -98,18 +106,18 @@ end @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.celldof_interior_indices(interpolation)]) - # Check for dirac delta property of interpolation - @testset "dirac delta property of dof $dof" for dof in 1:n_basefuncs - for k in 1:n_basefuncs - N_dof = shape_value(interpolation, coords[dof], k) - if k == dof - @test N_dof ≈ 1.0 - else - factor = interpolation isa Lagrange{RefQuadrilateral, 3} ? 200 : 4 - @test N_dof ≈ 0.0 atol = factor * eps(Float64) + # Check for dirac delta property of interpolation + @testset "dirac delta property of dof $dof" for dof in 1:n_basefuncs + for k in 1:n_basefuncs + N_dof = shape_value(interpolation, coords[dof], k) + if k == dof + @test N_dof ≈ 1.0 + else + factor = interpolation isa Lagrange{RefQuadrilateral, 3} ? 200 : 4 + @test N_dof ≈ 0.0 atol = factor * eps(Float64) + end end end - end # Test that facedof_indices(...) return in counter clockwise order (viewing from the outside) if interpolation isa Lagrange diff --git a/test/test_utils.jl b/test/test_utils.jl index d51fb8866c..c505d78acd 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -251,6 +251,28 @@ coords_on_faces(x, ::Serendipity{RefHexahedron, 2}) = check_equal_or_nan(a::Any, b::Any) = a==b || (isnan(a) && isnan(b)) check_equal_or_nan(a::Union{Tensor, Array}, b::Union{Tensor, Array}) = all(check_equal_or_nan.(a, b)) +# Hypercube is simply ⨂ᵈⁱᵐ Line :) +sample_random_point(::Type{Ferrite.RefHypercube{ref_dim}}) where {ref_dim} = Vec{ref_dim}(2.0 .* rand(Vec{ref_dim}) .- 1.0) +# Dirichlet type sampling +function sample_random_point(::Type{Ferrite.RefSimplex{ref_dim}}) where {ref_dim} + ξ = rand(ref_dim+1) + ξₜ = -log.(ξ) + return Vec{ref_dim}(ntuple(i->ξₜ[i], ref_dim) ./ sum(ξₜ)) +end +# Wedge = Triangle ⊗ Line +function sample_random_point(::Type{RefPrism}) + (ξ₁, ξ₂) = sample_random_point(RefTriangle) + ξ₃ = rand(Float64) + return Vec{3}((ξ₁, ξ₂, ξ₃)) +end +# TODO what to do here? The samplig is not uniform... +function sample_random_point(::Type{RefPyramid}) + ξ₃ = (1-1e-3)*rand(Float64) # Derivative is discontinuous at the top + # If we fix a z coordinate we get a Quad with extends (1-ξ₃) + (ξ₁, ξ₂) = (1.0 - ξ₃) .* rand(Vec{2}) + return Vec{3}((ξ₁, ξ₂, ξ₃)) +end + ###################################################### # Helpers for testing face_to_element_transformation # ###################################################### From b30c0bb7ccf5f5cbdb670dfe0bcf4f9bed99629e Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Mon, 6 Nov 2023 22:33:43 +0100 Subject: [PATCH 4/7] Make InterpolationInfo extendable by using outer constructor (#830) --- src/interpolations.jl | 89 ++++++++++++++++++------------------------- 1 file changed, 38 insertions(+), 51 deletions(-) diff --git a/src/interpolations.jl b/src/interpolations.jl index 4e1e70bf11..efe0f8b195 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -89,58 +89,44 @@ struct InterpolationInfo adjust_during_distribution::Bool n_copies::Int is_discontinuous::Bool - function InterpolationInfo(interpolation::InterpolationByDim{3}) - n_copies = 1 - if interpolation isa VectorizedInterpolation - n_copies = get_n_copies(interpolation) - interpolation = interpolation.ip - end - new( - [length(i) for i ∈ vertexdof_indices(interpolation)], - [length(i) for i ∈ edgedof_interior_indices(interpolation)], - [length(i) for i ∈ facedof_interior_indices(interpolation)], - length(celldof_interior_indices(interpolation)), - 3, - adjust_dofs_during_distribution(interpolation), - n_copies, - is_discontinuous(interpolation) - ) - end - function InterpolationInfo(interpolation::InterpolationByDim{2}) - n_copies = 1 - if interpolation isa VectorizedInterpolation - n_copies = get_n_copies(interpolation) - interpolation = interpolation.ip - end - new( - [length(i) for i ∈ vertexdof_indices(interpolation)], - Int[], - [length(i) for i ∈ facedof_interior_indices(interpolation)], - length(celldof_interior_indices(interpolation)), - 2, - adjust_dofs_during_distribution(interpolation), - n_copies, - is_discontinuous(interpolation) - ) - end - function InterpolationInfo(interpolation::InterpolationByDim{1}) - n_copies = 1 - if interpolation isa VectorizedInterpolation - n_copies = get_n_copies(interpolation) - interpolation = interpolation.ip - end - new( - [length(i) for i ∈ vertexdof_indices(interpolation)], - Int[], - Int[], - length(celldof_interior_indices(interpolation)), - 1, - adjust_dofs_during_distribution(interpolation), - n_copies, - is_discontinuous(interpolation) - ) - end end +function InterpolationInfo(interpolation::InterpolationByDim{3}, n_copies) + InterpolationInfo( + [length(i) for i ∈ vertexdof_indices(interpolation)], + [length(i) for i ∈ edgedof_interior_indices(interpolation)], + [length(i) for i ∈ facedof_interior_indices(interpolation)], + length(celldof_interior_indices(interpolation)), + 3, + adjust_dofs_during_distribution(interpolation), + n_copies, + is_discontinuous(interpolation) + ) +end +function InterpolationInfo(interpolation::InterpolationByDim{2}, n_copies) + InterpolationInfo( + [length(i) for i ∈ vertexdof_indices(interpolation)], + Int[], + [length(i) for i ∈ facedof_interior_indices(interpolation)], + length(celldof_interior_indices(interpolation)), + 2, + adjust_dofs_during_distribution(interpolation), + n_copies, + is_discontinuous(interpolation) + ) +end +function InterpolationInfo(interpolation::InterpolationByDim{1}, n_copies) + InterpolationInfo( + [length(i) for i ∈ vertexdof_indices(interpolation)], + Int[], + Int[], + length(celldof_interior_indices(interpolation)), + 1, + adjust_dofs_during_distribution(interpolation), + n_copies, + is_discontinuous(interpolation) + ) +end +InterpolationInfo(interpolation::Interpolation) = InterpolationInfo(interpolation, 1) # Some redundant information about the geometry of the reference cells. nfaces(::Interpolation{RefHypercube{dim}}) where {dim} = 2*dim @@ -1509,6 +1495,7 @@ end # Helper to get number of copies for DoF distribution get_n_copies(::VectorizedInterpolation{vdim}) where vdim = vdim +InterpolationInfo(ip::VectorizedInterpolation) = InterpolationInfo(ip.ip, get_n_copies(ip)) function getnbasefunctions(ipv::VectorizedInterpolation{vdim}) where vdim return vdim * getnbasefunctions(ipv.ip) From cbc0f70b464288575eca29a81cf4102731b54a7c Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Mon, 6 Nov 2023 22:34:14 +0100 Subject: [PATCH 5/7] More bugfixes and regression tests for the topology. (#758) --- src/Grid/topology.jl | 18 ++++++++---------- test/test_grid_dofhandler_vtk.jl | 27 +++++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 10 deletions(-) diff --git a/src/Grid/topology.jl b/src/Grid/topology.jl index ed5f1e26ac..a34b3acadb 100644 --- a/src/Grid/topology.jl +++ b/src/Grid/topology.jl @@ -260,17 +260,15 @@ function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, vertexidx:: cellid, local_vertexid = vertexidx[1], vertexidx[2] cell_vertices = vertices(getcells(grid,cellid)) global_vertexid = cell_vertices[local_vertexid] - if include_self - vertex_to_cell = top.vertex_to_cell[global_vertexid] - self_reference_local = Vector{VertexIndex}(undef,length(vertex_to_cell)) - for (i,cellid) in enumerate(vertex_to_cell) - local_vertex = VertexIndex(cellid,findfirst(x->x==global_vertexid,vertices(getcells(grid,cellid)))::Int) - self_reference_local[i] = local_vertex - end - return [top.vertex_vertex_neighbor[global_vertexid].neighbor_info; self_reference_local] - else - return top.vertex_vertex_neighbor[global_vertexid].neighbor_info + vertex_to_cell = top.vertex_to_cell[global_vertexid] + self_reference_local = Vector{VertexIndex}() + sizehint!(self_reference_local, length(vertex_to_cell)) + for (i,cellid) in enumerate(vertex_to_cell) + local_vertex = VertexIndex(cellid,findfirst(x->x==global_vertexid,vertices(getcells(grid,cellid)))::Int) + !include_self && local_vertex == vertexidx && continue + push!(self_reference_local, local_vertex) end + return self_reference_local end function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid{3}, edgeidx::EdgeIndex, include_self=false) diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 05ca612075..f29ad7969f 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/test/test_grid_dofhandler_vtk.jl @@ -286,9 +286,13 @@ end linegrid = generate_grid(Line,(3,)) linetopo = ExclusiveTopology(linegrid) @test linetopo.vertex_vertex_neighbor[1,2] == Ferrite.EntityNeighborhood(VertexIndex(2,1)) + @test getneighborhood(linetopo, linegrid, VertexIndex(1,2)) == [VertexIndex(2,1)] @test linetopo.vertex_vertex_neighbor[2,1] == Ferrite.EntityNeighborhood(VertexIndex(1,2)) + @test getneighborhood(linetopo, linegrid, VertexIndex(2,1)) == [VertexIndex(1,2)] @test linetopo.vertex_vertex_neighbor[2,2] == Ferrite.EntityNeighborhood(VertexIndex(3,1)) + @test getneighborhood(linetopo, linegrid, VertexIndex(2,2)) == [VertexIndex(3,1)] @test linetopo.vertex_vertex_neighbor[3,1] == Ferrite.EntityNeighborhood(VertexIndex(2,2)) + @test getneighborhood(linetopo, linegrid, VertexIndex(3,1)) == [VertexIndex(2,2)] linefaceskeleton = Ferrite.faceskeleton(linetopo, linegrid) @test length(linefaceskeleton) == 4 @@ -306,31 +310,52 @@ end #test vertex neighbors maps cellid and local vertex id to neighbor id and neighbor local vertex id @test topology.vertex_vertex_neighbor[1,3] == Ferrite.EntityNeighborhood(VertexIndex(4,1)) @test topology.vertex_vertex_neighbor[2,4] == Ferrite.EntityNeighborhood(VertexIndex(3,2)) + @test Set(getneighborhood(topology, quadgrid, VertexIndex(2,4))) == Set([VertexIndex(1,3), VertexIndex(3,2), VertexIndex(4,1)]) @test topology.vertex_vertex_neighbor[3,3] == Ferrite.EntityNeighborhood(VertexIndex(6,1)) @test topology.vertex_vertex_neighbor[3,2] == Ferrite.EntityNeighborhood(VertexIndex(2,4)) @test topology.vertex_vertex_neighbor[4,1] == Ferrite.EntityNeighborhood(VertexIndex(1,3)) @test topology.vertex_vertex_neighbor[4,4] == Ferrite.EntityNeighborhood(VertexIndex(5,2)) @test topology.vertex_vertex_neighbor[5,2] == Ferrite.EntityNeighborhood(VertexIndex(4,4)) @test topology.vertex_vertex_neighbor[6,1] == Ferrite.EntityNeighborhood(VertexIndex(3,3)) + @test isempty(getneighborhood(topology, quadgrid, VertexIndex(2,2))) + @test length(getneighborhood(topology, quadgrid, VertexIndex(2,4))) == 3 #test face neighbor maps cellid and local face id to neighbor id and neighbor local face id @test topology.face_face_neighbor[1,2] == Ferrite.EntityNeighborhood(FaceIndex(2,4)) + @test getneighborhood(topology, quadgrid, FaceIndex(1,2)) == [FaceIndex(2,4)] @test topology.face_face_neighbor[1,3] == Ferrite.EntityNeighborhood(FaceIndex(3,1)) + @test getneighborhood(topology, quadgrid, FaceIndex(1,3)) == [FaceIndex(3,1)] @test topology.face_face_neighbor[2,3] == Ferrite.EntityNeighborhood(FaceIndex(4,1)) + @test getneighborhood(topology, quadgrid, FaceIndex(2,3)) == [FaceIndex(4,1)] @test topology.face_face_neighbor[2,4] == Ferrite.EntityNeighborhood(FaceIndex(1,2)) + @test getneighborhood(topology, quadgrid, FaceIndex(2,4)) == [FaceIndex(1,2)] @test topology.face_face_neighbor[3,1] == Ferrite.EntityNeighborhood(FaceIndex(1,3)) + @test getneighborhood(topology, quadgrid, FaceIndex(3,1)) == [FaceIndex(1,3)] @test topology.face_face_neighbor[3,2] == Ferrite.EntityNeighborhood(FaceIndex(4,4)) + @test getneighborhood(topology, quadgrid, FaceIndex(3,2)) == [FaceIndex(4,4)] @test topology.face_face_neighbor[3,3] == Ferrite.EntityNeighborhood(FaceIndex(5,1)) + @test getneighborhood(topology, quadgrid, FaceIndex(3,3)) == [FaceIndex(5,1)] @test topology.face_face_neighbor[4,1] == Ferrite.EntityNeighborhood(FaceIndex(2,3)) + @test getneighborhood(topology, quadgrid, FaceIndex(4,1)) == [FaceIndex(2,3)] @test topology.face_face_neighbor[4,3] == Ferrite.EntityNeighborhood(FaceIndex(6,1)) + @test getneighborhood(topology, quadgrid, FaceIndex(4,3)) == [FaceIndex(6,1)] @test topology.face_face_neighbor[4,4] == Ferrite.EntityNeighborhood(FaceIndex(3,2)) + @test getneighborhood(topology, quadgrid, FaceIndex(1,2)) == [FaceIndex(2,4)] @test topology.face_face_neighbor[5,1] == Ferrite.EntityNeighborhood(FaceIndex(3,3)) + @test getneighborhood(topology, quadgrid, FaceIndex(5,1)) == [FaceIndex(3,3)] @test topology.face_face_neighbor[5,2] == Ferrite.EntityNeighborhood(FaceIndex(6,4)) + @test getneighborhood(topology, quadgrid, FaceIndex(5,2)) == [FaceIndex(6,4)] @test topology.face_face_neighbor[5,3] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[]) + @test getneighborhood(topology, quadgrid, FaceIndex(5,3)) == FaceIndex[] @test topology.face_face_neighbor[5,4] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[]) + @test getneighborhood(topology, quadgrid, FaceIndex(5,4)) == FaceIndex[] @test topology.face_face_neighbor[6,1] == Ferrite.EntityNeighborhood(FaceIndex(4,3)) + @test getneighborhood(topology, quadgrid, FaceIndex(6,1)) == [FaceIndex(4,3)] @test topology.face_face_neighbor[6,2] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[]) + @test getneighborhood(topology, quadgrid, FaceIndex(6,2)) == FaceIndex[] @test topology.face_face_neighbor[6,3] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[]) + @test getneighborhood(topology, quadgrid, FaceIndex(6,3)) == FaceIndex[] @test topology.face_face_neighbor[6,4] == Ferrite.EntityNeighborhood(FaceIndex(5,2)) + @test getneighborhood(topology, quadgrid, FaceIndex(6,4)) == [FaceIndex(5,2)] # (8) # (7) +-----+-----+(9) # | 3 | 4 | @@ -349,8 +374,10 @@ end topology = ExclusiveTopology(hexgrid) @test topology.edge_edge_neighbor[1,11] == Ferrite.EntityNeighborhood(EdgeIndex(4,9)) @test Set(getneighborhood(topology,hexgrid,EdgeIndex(1,11),true)) == Set([EdgeIndex(4,9),EdgeIndex(2,12),EdgeIndex(3,10),EdgeIndex(1,11)]) + @test Set(getneighborhood(topology,hexgrid,EdgeIndex(1,11),false)) == Set([EdgeIndex(4,9),EdgeIndex(2,12),EdgeIndex(3,10)]) @test topology.edge_edge_neighbor[2,12] == Ferrite.EntityNeighborhood(EdgeIndex(3,10)) @test Set(getneighborhood(topology,hexgrid,EdgeIndex(2,12),true)) == Set([EdgeIndex(3,10),EdgeIndex(1,11),EdgeIndex(4,9),EdgeIndex(2,12)]) + @test Set(getneighborhood(topology,hexgrid,EdgeIndex(2,12),false)) == Set([EdgeIndex(3,10),EdgeIndex(1,11),EdgeIndex(4,9)]) @test topology.edge_edge_neighbor[3,10] == Ferrite.EntityNeighborhood(EdgeIndex(2,12)) @test topology.edge_edge_neighbor[4,9] == Ferrite.EntityNeighborhood(EdgeIndex(1,11)) @test getneighborhood(topology,hexgrid,FaceIndex((1,3))) == [FaceIndex((2,5))] From 13f2aee53d339f368d33b7acbf9ae72f0960b77f Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Mon, 6 Nov 2023 23:54:57 +0100 Subject: [PATCH 6/7] Helper script for convergence rates for Poisson problems (#640) --- src/exports.jl | 1 + src/interpolations.jl | 1 + .../test_simple_scalar_convergence.jl | 216 ++++++++++++++++++ test/runtests.jl | 5 + 4 files changed, 223 insertions(+) create mode 100644 test/integration/test_simple_scalar_convergence.jl diff --git a/src/exports.jl b/src/exports.jl index 09a0ab1e5e..6ab0507a7a 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -17,6 +17,7 @@ export DiscontinuousLagrange, Serendipity, getnbasefunctions, + getrefshape, # Quadrature QuadratureRule, diff --git a/src/interpolations.jl b/src/interpolations.jl index efe0f8b195..8aa03dd533 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -1268,6 +1268,7 @@ end ####################################### # Taken from https://defelement.com/elements/bubble-enriched-lagrange.html getnbasefunctions(::BubbleEnrichedLagrange{RefTriangle,1}) = 4 +adjust_dofs_during_distribution(::BubbleEnrichedLagrange{RefTriangle,1}) = false vertexdof_indices(::BubbleEnrichedLagrange{RefTriangle,1}) = ((1,), (2,), (3,)) facedof_indices(::BubbleEnrichedLagrange{RefTriangle,1}) = ((1,2), (2,3), (3,1)) diff --git a/test/integration/test_simple_scalar_convergence.jl b/test/integration/test_simple_scalar_convergence.jl new file mode 100644 index 0000000000..d6387fc0d3 --- /dev/null +++ b/test/integration/test_simple_scalar_convergence.jl @@ -0,0 +1,216 @@ +using Ferrite, Test +import Ferrite: getdim, default_interpolation + +module ConvergenceTestHelper + +using Ferrite, SparseArrays, ForwardDiff, Test +import LinearAlgebra: diag + +get_geometry(::Ferrite.Interpolation{RefLine}) = Line +get_geometry(::Ferrite.Interpolation{RefQuadrilateral}) = Quadrilateral +get_geometry(::Ferrite.Interpolation{RefTriangle}) = Triangle +get_geometry(::Ferrite.Interpolation{RefPrism}) = Wedge +get_geometry(::Ferrite.Interpolation{RefHexahedron}) = Hexahedron +get_geometry(::Ferrite.Interpolation{RefTetrahedron}) = Tetrahedron +get_geometry(::Ferrite.Interpolation{RefPyramid}) = Pyramid + +get_quadrature_order(::Lagrange{shape, order}) where {shape, order} = 2*order +get_quadrature_order(::Serendipity{shape, order}) where {shape, order} = 2*order +get_quadrature_order(::CrouzeixRaviart{shape, order}) where {shape, order} = 2*order+1 +get_quadrature_order(::BubbleEnrichedLagrange{shape, order}) where {shape, order} = 2*order + +get_num_elements(::Ferrite.Interpolation{shape, 1}) where {shape} = 21 +get_num_elements(::Ferrite.Interpolation{shape, 2}) where {shape} = 7 +get_num_elements(::Ferrite.Interpolation{RefHexahedron, 1}) = 11 +get_num_elements(::Ferrite.Interpolation{RefHexahedron, 2}) = 4 +get_num_elements(::Ferrite.Interpolation{shape, 3}) where {shape} = 8 +get_num_elements(::Ferrite.Interpolation{shape, 4}) where {shape} = 5 +get_num_elements(::Ferrite.Interpolation{shape, 5}) where {shape} = 3 + +analytical_solution(x) = prod(cos, x*π/2) +analytical_rhs(x) = -Tensors.laplace(analytical_solution,x) + +# Standard assembly copy pasta for Poisson problem +function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues, coords) + n_basefuncs = getnbasefunctions(cellvalues) + ## Reset to 0 + fill!(Ke, 0) + fill!(fe, 0) + ## Loop over quadrature points + for q_point in 1:getnquadpoints(cellvalues) + ## Get the quadrature weight + dΩ = getdetJdV(cellvalues, q_point) + x = spatial_coordinate(cellvalues, q_point, coords) + ## Loop over test shape functions + for i in 1:n_basefuncs + δu = shape_value(cellvalues, q_point, i) + ∇δu = shape_gradient(cellvalues, q_point, i) + ## Add contribution to fe + fe[i] += analytical_rhs(x) * δu * dΩ + ## Loop over trial shape functions + for j in 1:n_basefuncs + ∇u = shape_gradient(cellvalues, q_point, j) + ## Add contribution to Ke + Ke[i, j] += (∇δu ⋅ ∇u) * dΩ + end + end + end + return Ke, fe +end + +# Standard assembly copy pasta for Poisson problem +function assemble_global(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler) + ## Allocate the element stiffness matrix and element force vector + n_basefuncs = getnbasefunctions(cellvalues) + Ke = zeros(n_basefuncs, n_basefuncs) + fe = zeros(n_basefuncs) + ## Allocate global force vector f + f = zeros(ndofs(dh)) + ## Create an assembler + assembler = start_assemble(K, f) + ## Loop over all cels + for cell in CellIterator(dh) + ## Reinitialize cellvalues for this cell + reinit!(cellvalues, cell) + coords = getcoordinates(cell) + ## Compute element contribution + assemble_element!(Ke, fe, cellvalues, coords) + ## Assemble Ke and fe into K and f + assemble!(assembler, celldofs(cell), Ke, fe) + end + return K, f +end + +# Compute norms +function check_and_compute_convergence_norms(dh, u, cellvalues, testatol) + L2norm = 0.0 + ∇L2norm = 0.0 + L∞norm = 0.0 + for cell in CellIterator(dh) + reinit!(cellvalues, cell) + n_basefuncs = getnbasefunctions(cellvalues) + coords = getcoordinates(cell) + uₑ = u[celldofs(cell)] + for q_point in 1:getnquadpoints(cellvalues) + dΩ = getdetJdV(cellvalues, q_point) + x = spatial_coordinate(cellvalues, q_point, coords) + uₐₙₐ = prod(cos, x*π/2) + uₐₚₚᵣₒₓ = function_value(cellvalues, q_point, uₑ) + L∞norm = max(L∞norm, norm(uₐₙₐ-uₐₚₚᵣₒₓ)) + L2norm += norm(uₐₙₐ-uₐₚₚᵣₒₓ)^2*dΩ + + ∇uₐₙₐ = gradient(x-> prod(cos, x*π/2), x) + ∇uₐₚₚᵣₒₓ = function_gradient(cellvalues, q_point, uₑ) + ∇L2norm += norm(∇uₐₙₐ-∇uₐₚₚᵣₒₓ)^2*dΩ + + # Pointwise convergence + @test uₐₙₐ ≈ uₐₚₚᵣₒₓ atol=testatol + end + end + √(L2norm), √(∇L2norm), L∞norm +end + +# Assemble and solve +function solve(dh, ch, cellvalues) + K, f = assemble_global(cellvalues, create_sparsity_pattern(dh), dh); + apply!(K, f, ch) + u = K \ f; +end + +function setup_poisson_problem(grid, interpolation, interpolation_geo, qr) + # Construct Ferrite stuff + dh = DofHandler(grid) + add!(dh, :u, interpolation) + close!(dh); + + ch = ConstraintHandler(dh); + ∂Ω = union( + values(grid.facesets)... + ); + dbc = Dirichlet(:u, ∂Ω, (x, t) -> analytical_solution(x)) + add!(ch, dbc); + close!(ch); + + cellvalues = CellValues(qr, interpolation, interpolation_geo); + + return dh, ch, cellvalues +end + +end # module ConvergenceTestHelper + +# These test only for convergence within margins +@testset "convergence analysis" begin + @testset "$interpolation" for interpolation in ( + Lagrange{RefTriangle, 3}(), + Lagrange{RefTriangle, 4}(), + Lagrange{RefTriangle, 5}(), + Lagrange{RefHexahedron, 1}(), + Lagrange{RefTetrahedron, 1}(), + Lagrange{RefPrism, 1}(), + Lagrange{RefPyramid, 1}(), + # + Serendipity{RefQuadrilateral, 2}(), + Serendipity{RefHexahedron, 2}(), + # + BubbleEnrichedLagrange{RefTriangle, 1}(), + # + CrouzeixRaviart{RefTriangle, 1}(), + ) + # Generate a grid ... + geometry = ConvergenceTestHelper.get_geometry(interpolation) + interpolation_geo = default_interpolation(geometry) + N = ConvergenceTestHelper.get_num_elements(interpolation) + grid = generate_grid(geometry, ntuple(x->N, getdim(geometry))); + # ... a suitable quadrature rule ... + qr_order = ConvergenceTestHelper.get_quadrature_order(interpolation) + qr = QuadratureRule{getrefshape(interpolation)}(qr_order) + # ... and then pray to the gods of convergence. + dh, ch, cellvalues = ConvergenceTestHelper.setup_poisson_problem(grid, interpolation, interpolation_geo, qr) + u = ConvergenceTestHelper.solve(dh, ch, cellvalues) + ConvergenceTestHelper.check_and_compute_convergence_norms(dh, u, cellvalues, 1e-2) + end +end + +# These test also for correct convergence rates +@testset "convergence rate" begin + @testset "$interpolation" for interpolation in ( + Lagrange{RefLine, 1}(), + Lagrange{RefLine, 2}(), + Lagrange{RefQuadrilateral, 1}(), + Lagrange{RefQuadrilateral, 2}(), + Lagrange{RefQuadrilateral, 3}(), + Lagrange{RefTriangle, 1}(), + Lagrange{RefTriangle, 2}(), + Lagrange{RefHexahedron, 2}(), + Lagrange{RefTetrahedron, 2}(), + Lagrange{RefPrism, 2}(), + ) + # Generate a grid ... + geometry = ConvergenceTestHelper.get_geometry(interpolation) + interpolation_geo = default_interpolation(geometry) + # "Coarse case" + N₁ = ConvergenceTestHelper.get_num_elements(interpolation) + grid = generate_grid(geometry, ntuple(x->N₁, getdim(geometry))); + # ... a suitable quadrature rule ... + qr_order = ConvergenceTestHelper.get_quadrature_order(interpolation) + qr = QuadratureRule{getrefshape(interpolation)}(qr_order) + # ... and then pray to the gods of convergence. + dh, ch, cellvalues = ConvergenceTestHelper.setup_poisson_problem(grid, interpolation, interpolation_geo, qr) + u = ConvergenceTestHelper.solve(dh, ch, cellvalues) + L2₁, H1₁, _ = ConvergenceTestHelper.check_and_compute_convergence_norms(dh, u, cellvalues, 1e-2) + + # "Fine case" + N₂ = 2*N₁ + grid = generate_grid(geometry, ntuple(x->N₂, getdim(geometry))); + # ... a suitable quadrature rule ... + qr_order = ConvergenceTestHelper.get_quadrature_order(interpolation) + qr = QuadratureRule{getrefshape(interpolation)}(qr_order) + # ... and then pray to the gods of convergence. + dh, ch, cellvalues = ConvergenceTestHelper.setup_poisson_problem(grid, interpolation, interpolation_geo, qr) + u = ConvergenceTestHelper.solve(dh, ch, cellvalues) + L2₂, H1₂, _ = ConvergenceTestHelper.check_and_compute_convergence_norms(dh, u, cellvalues, 5e-3) + + @test -(log(L2₂)-log(L2₁))/(log(N₂)-log(N₁)) ≈ Ferrite.getorder(interpolation)+1 atol=0.1 + @test -(log(H1₂)-log(H1₁))/(log(N₂)-log(N₁)) ≈ Ferrite.getorder(interpolation) atol=0.1 + end +end diff --git a/test/runtests.jl b/test/runtests.jl index d11ed6d582..fe3c1692c9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,8 @@ if HAS_EXTENSIONS && MODULE_CAN_BE_TYPE_PARAMETER end include("test_utils.jl") + +# Unit tests include("test_interpolations.jl") include("test_cellvalues.jl") include("test_facevalues.jl") @@ -38,3 +40,6 @@ include("test_deprecations.jl") HAS_EXTENSIONS && include("blockarrays.jl") include("test_examples.jl") @test all(x -> isdefined(Ferrite, x), names(Ferrite)) # Test that all exported symbols are defined + +# Integration tests +include("integration/test_simple_scalar_convergence.jl") From d0f0e07456ba3f82b34f1fba7ecaa8c0dd41207d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20K=C3=B6hler?= Date: Tue, 7 Nov 2023 01:10:04 +0100 Subject: [PATCH 7/7] add isequal and hash methods to check tuple in BoundaryIndex Sets without type piracy (#835) --- CHANGELOG.md | 5 +++++ src/Grid/grid.jl | 7 +++++-- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b3d483c405..6c1af9946c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -319,6 +319,10 @@ more discussion). **To upgrade** replace the quadrature rule passed to `FaceValues` with a `FaceQuadratureRule`. +- Checking if a face `(ele_id, local_face_id) ∈ faceset` has been previously implemented + by type piracy. In order to be invariant to the underlying `Set` datatype as well as + omitting type piracy, ([#835][github-835]) implemented `isequal` and `hash` for `BoundaryIndex` datatypes. + ### Deprecated - The rarely (if ever) used methods of `function_value`, `function_gradient`, @@ -867,3 +871,4 @@ poking into Ferrite internals: [github-749]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/749 [github-753]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/753 [github-756]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/756 +[github-835]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/835 diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index d3b7bc68cd..cab1f7409b 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -765,8 +765,11 @@ for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex) #To be able to do a,b = faceidx Base.iterate(I::($INDEX), state::Int=1) = (state==3) ? nothing : (I[state], state+1) - #For (cellid, faceidx) in faceset - Base.in(v::Tuple{Int, Int}, s::Set{$INDEX}) = in($INDEX(v), s) + # Necessary to check if, e.g. `(cellid, faceidx) in faceset` + Base.isequal(x::$INDEX, y::$INDEX) = x.idx == y.idx + Base.isequal(x::Tuple{Int, Int}, y::$INDEX) = x[1] == y.idx[1] && x[2] == y.idx[2] + Base.isequal(y::$INDEX, x::Tuple{Int, Int}) = x[1] == y.idx[1] && x[2] == y.idx[2] + Base.hash(x::$INDEX, h::UInt) = hash(x.idx, h) end end