From 7c6e055d84d911b6e851121773b38c6a2ba585c7 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Thu, 6 Oct 2022 11:15:17 +0200 Subject: [PATCH] Default to constraining all components for (Periodic)Dirichlet. (#509) This patch changes the default components to prescribe for (Periodic)Dirichlet from component one to all components. This change does not affect scalar problems, and judging by tests, does not really affect vector valued problems either since it was already rather strange to rely on the default behavior of prescribing only the first component. Fixes #506. --- CHANGELOG.md | 6 ++ src/Dofs/ConstraintHandler.jl | 109 ++++++++++++++++++++-------------- src/Dofs/DofHandler.jl | 2 +- src/Dofs/MixedDofHandler.jl | 6 +- test/test_constraints.jl | 75 +++++++++++++++++++++++ 5 files changed, 147 insertions(+), 51 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index eb42d83749..bfb1224fd9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Changed + - The default components to constrain in `Dirichlet` and `PeriodicDirichlet` have changed + from component 1 to all components of the field. For scalar problems this has no effect. + ([#506][github-506], [#509][github-509]) ## [0.3.8] - 2022-10-05 ### Added @@ -123,6 +127,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 [github-501]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/501 [github-503]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/503 [github-505]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/505 +[github-506]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/506 +[github-509]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/509 [Unreleased]: https://github.com/Ferrite-FEM/Ferrite.jl/compare/v0.3.8...HEAD [0.3.8]: https://github.com/Ferrite-FEM/Ferrite.jl/compare/v0.3.7...v0.3.8 diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 9a2cebfbc0..bf202098fc 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -1,26 +1,31 @@ # abstract type Constraint end """ - Dirichlet(u, ∂Ω, f) - Dirichlet(u, ∂Ω, f, component) + Dirichlet(u::Symbol, ∂Ω::Set, f::Function, components=nothing) Create a Dirichlet boundary condition on `u` on the `∂Ω` part of the boundary. `f` is a function that takes two arguments, `x` and `t` where `x` is the spatial coordinate and `t` is the current time, -and returns the prescribed value. For example, here we create a +and returns the prescribed value. `components` specify the components +of `u` that are prescribed by this condition. By default all components +of `u` are prescribed. + +For example, here we create a Dirichlet condition for the `:u` field, on the faceset called `∂Ω` and the value given by the `sin` function: -```julia -dbc = Dirichlet(:u, ∂Ω, (x, t) -> sin(t)) -``` +*Examples* +```jldoctest +# Obtain the faceset from the grid +∂Ω = getfaceset(grid, "boundary-1") -If `:u` is a vector field we can specify which component the condition -should be applied to by specifying `component`. `component` can be given -either as an integer, or as a vector, for example: +# Prescribe scalar field :s on ∂Ω to sin(t) +dbc = Dirichlet(:s, ∂Ω, (x, t) -> sin(t)) -```julia -dbc = Dirichlet(:u, ∂Ω, (x, t) -> sin(t), 1) # applied to component 1 -dbc = Dirichlet(:u, ∂Ω, (x, t) -> sin(t), [1, 3]) # applied to component 1 and 3 +# Prescribe all components of vector field :v on ∂Ω to 0 +dbc = Dirichlet(:v, ∂Ω, (x, t) -> 0 * x) + +# Prescribe component 2 and 3 of vector field :v on ∂Ω to [sin(t), cos(t)] +dbc = Dirichlet(:v, ∂Ω, (x, t) -> [sin(t), cos(t)], [2, 3]) ``` `Dirichlet` boundary conditions are added to a [`ConstraintHandler`](@ref) @@ -34,14 +39,18 @@ struct Dirichlet # <: Constraint local_face_dofs::Vector{Int} local_face_dofs_offset::Vector{Int} end -function Dirichlet(field_name::Symbol, faces::Set, f::Function, components=1) +function Dirichlet(field_name::Symbol, faces::Set, f::Function, components=nothing) return Dirichlet(f, copy(faces), field_name, __to_components(components), Int[], Int[]) end +# components=nothing is default and means that all components should be constrained +# but since number of components isn't known here it will be populated in add! +__to_components(::Nothing) = Int[] function __to_components(c) components = convert(Vector{Int}, vec(collect(Int, c))) - issorted(components) || error("components not sorted: $components") - allunique(components) || error("components not unique: $components") + isempty(components) && error("components are empty: $c") + issorted(components) || error("components not sorted: $c") + allunique(components) || error("components not unique: $c") return components end @@ -191,33 +200,29 @@ function close!(ch::ConstraintHandler) return ch end -function dbc_check(ch::ConstraintHandler, dbc::Dirichlet) - # check input - dbc.field_name in getfieldnames(ch.dh) || throw(ArgumentError("field $(dbc.field_name) does not exist in DofHandler, existing fields are $(getfieldnames(ch.dh))")) - #TODO FIX!! - #for component in dbc.components - # 0 < component <= ndim(ch.dh, dbc.field_name) || error("component $component is not within the range of field $field which has $(ndim(ch.dh, field)) dimensions") - #end - if length(dbc.faces) == 0 - @warn("added Dirichlet Boundary Condition to set containing 0 entities") - end -end - """ add!(ch::ConstraintHandler, dbc::Dirichlet) Add a `Dirichlet` boundary condition to the `ConstraintHandler`. """ function add!(ch::ConstraintHandler, dbc::Dirichlet) - dbc_check(ch, dbc) + if length(dbc.faces) == 0 + @warn("adding Dirichlet Boundary Condition to set containing 0 entities") + end celltype = getcelltype(ch.dh.grid) @assert isconcretetype(celltype) - field_idx = find_field(ch.dh, dbc.field_name) # Extract stuff for the field - interpolation = getfieldinterpolation(ch.dh, field_idx)#ch.dh.field_interpolations[field_idx] - field_dim = getfielddim(ch.dh, field_idx)#ch.dh.field_dims[field_idx] # TODO: I think we don't need to extract these here ... + field_idx = find_field(ch.dh, dbc.field_name) # throws if name not found + interpolation = getfieldinterpolation(ch.dh, field_idx) + field_dim = getfielddim(ch.dh, field_idx) + + if !all(c -> 0 < c <= field_dim, dbc.components) + error("components $(dbc.components) not within range of field :$(dbc.field_name) ($(field_dim) dimension(s))") + end + # Empty components means constrain them all + isempty(dbc.components) && append!(dbc.components, 1:field_dim) if eltype(dbc.faces)==Int #Special case when dbc.faces is a nodeset bcvalue = BCValues(interpolation, default_interpolation(celltype), FaceIndex) #Not used by node bcs, but still have to pass it as an argument @@ -816,11 +821,18 @@ function add!(ch::ConstraintHandler, fh::FieldHandler, dbc::Dirichlet) celltype = getcelltype(ch.dh.grid, first(fh.cellset)) #Assume same celltype of all cells in fh.cellset - field_idx = find_field(fh, dbc.field_name) # Extract stuff for the field + field_idx = find_field(fh, dbc.field_name) interpolation = getfieldinterpolations(fh)[field_idx] field_dim = getfielddims(fh)[field_idx] + if !all(c -> 0 < c <= field_dim, dbc.components) + error("components $(dbc.components) not within range of field :$(dbc.field_name) ($(field_dim) dimension(s))") + end + + # Empty components means constrain them all + isempty(dbc.components) && append!(dbc.components, 1:field_dim) + if eltype(dbc.faces)==Int #Special case when dbc.faces is a nodeset bcvalue = BCValues(interpolation, default_interpolation(celltype), FaceIndex) #Not used by node bcs, but still have to pass it as an argument else @@ -868,15 +880,15 @@ struct PeriodicFacePair end """ - PeriodicDirichlet(u, face_mapping, component=1) - PeriodicDirichlet(u, face_mapping, R::AbstractMatrix, component=1) - PeriodicDirichlet(u, face_mapping, f::Function, component=1) + PeriodicDirichlet(u::Symbol, face_mapping, components=nothing) + PeriodicDirichlet(u::Symbol, face_mapping, R::AbstractMatrix, components=nothing) + PeriodicDirichlet(u::Symbol, face_mapping, f::Function, components=nothing) Create a periodic Dirichlet boundary condition for the field `u` on the face-pairs given in `face_mapping`. The mapping can be computed with [`collect_periodic_faces`](@ref). The constraint ensures that degrees-of-freedom on the mirror face are constrained to the -corresponding degrees-of-freedom on the image face. The condition is imposed in a -strong sense, and requires (i) a periodic domain and (ii) a periodic mesh. +corresponding degrees-of-freedom on the image face. `components` specify the components of +`u` that are prescribed by this condition. By default all components of `u` are prescribed. If the mapping is not aligned with the coordinate axis (e.g. rotated) a rotation matrix `R` should be passed to the constructor. This matrix rotates dofs on the mirror face to the @@ -898,25 +910,20 @@ struct PeriodicDirichlet end # Default to no inhomogeneity function/rotation -PeriodicDirichlet(fn::Symbol, fp::Union{Vector{<:Pair},Vector{PeriodicFacePair}}, c=1) = +PeriodicDirichlet(fn::Symbol, fp::Union{Vector{<:Pair},Vector{PeriodicFacePair}}, c=nothing) = PeriodicDirichlet(fn, fp, nothing, c) # Basic constructor for the simple case where face_map will be populated in # add!(::ConstraintHandler, ...) instead -function PeriodicDirichlet(fn::Symbol, fp::Vector{<:Pair}, f::Union{Function,Nothing}, c=1) +function PeriodicDirichlet(fn::Symbol, fp::Vector{<:Pair}, f::Union{Function,Nothing}, c=nothing) face_map = PeriodicFacePair[] # This will be populated in add!(::ConstraintHandler, ...) instead return PeriodicDirichlet(fn, __to_components(c), fp, face_map, f, nothing) end -function PeriodicDirichlet(fn::Symbol, fm::Vector{PeriodicFacePair}, f_or_r::Union{AbstractMatrix,Function,Nothing}, c=1) +function PeriodicDirichlet(fn::Symbol, fm::Vector{PeriodicFacePair}, f_or_r::Union{AbstractMatrix,Function,Nothing}, c=nothing) f = f_or_r isa Function ? f_or_r : nothing rotation_matrix = f_or_r isa AbstractMatrix ? f_or_r : nothing components = __to_components(c) - if rotation_matrix !== nothing - if !(length(components) == size(rotation_matrix, 1) == size(rotation_matrix, 2)) - error("size of rotation matrix does not match the number of components") - end - end return PeriodicDirichlet(fn, components, Pair{String,String}[], fm, f, rotation_matrix) end @@ -931,6 +938,14 @@ function add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet) field_idx = find_field(ch.dh, pdbc.field_name) interpolation = getfieldinterpolation(ch.dh, field_idx) field_dim = getfielddim(ch.dh, field_idx) + + if !all(c -> 0 < c <= field_dim, pdbc.components) + error("components $(pdbc.components) not within range of field :$(pdbc.field_name) ($(field_dim) dimension(s))") + end + + # Empty components means constrain them all + isempty(pdbc.components) && append!(pdbc.components, 1:field_dim) + if pdbc.rotation_matrix === nothing dof_map_t = Int iterator_f = identity @@ -940,7 +955,9 @@ function add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet) error("legacy mode not supported with rotations") end nc = length(pdbc.components) - @assert nc == size(pdbc.rotation_matrix, 1) == size(pdbc.rotation_matrix, 2) # Verified in constructor + if !(nc == size(pdbc.rotation_matrix, 1) == size(pdbc.rotation_matrix, 2)) + error("size of rotation matrix does not match the number of components") + end if nc !== field_dim error("rotations currently only supported when all components are periodic") end diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index 1483a29b12..83cd9bebac 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -55,7 +55,7 @@ getfieldnames(dh::AbstractDofHandler) = dh.field_names ndim(dh::AbstractDofHandler, field_name::Symbol) = dh.field_dims[find_field(dh, field_name)] function find_field(dh::DofHandler, field_name::Symbol) j = findfirst(i->i == field_name, dh.field_names) - j == 0 && error("did not find field $field_name") + j === nothing && error("could not find field :$field_name in DofHandler (existing fields: $(getfieldnames(dh)))") return j end diff --git a/src/Dofs/MixedDofHandler.jl b/src/Dofs/MixedDofHandler.jl index 6b02a01383..d4cf134aea 100644 --- a/src/Dofs/MixedDofHandler.jl +++ b/src/Dofs/MixedDofHandler.jl @@ -102,9 +102,7 @@ Returns the union of all the fields. Can be used as an iterable over all the fie function getfieldnames(dh::MixedDofHandler) fieldnames = Vector{Symbol}() for fh in dh.fieldhandlers - for name in getfieldnames(fh) - push!(fieldnames, name) - end + append!(fieldnames, getfieldnames(fh)) end return unique!(fieldnames) end @@ -441,7 +439,7 @@ create_symmetric_sparsity_pattern(dh::MixedDofHandler) = Symmetric(_create_spars function find_field(fh::FieldHandler, field_name::Symbol) j = findfirst(i->i == field_name, getfieldnames(fh)) - j === nothing && error("did not find field $field_name") + j === nothing && error("could not find field :$field_name in FieldHandler (existing fields: $(getfieldnames(fh)))") return j end diff --git a/test/test_constraints.jl b/test/test_constraints.jl index 74ac57405b..0fc7ade75c 100644 --- a/test/test_constraints.jl +++ b/test/test_constraints.jl @@ -1,4 +1,79 @@ # misc constraint tests + +@testset "constructors and error checking" begin + grid = generate_grid(Triangle, (2, 2)) + Γ = getfaceset(grid, "left") + face_map = collect_periodic_faces(grid, "left", "right") + dh = DofHandler(grid) + push!(dh, :s, 1) + push!(dh, :v, 2) + close!(dh) + ch = ConstraintHandler(dh) + + # Dirichlet + @test_throws ErrorException("components are empty: $(Int)[]") Dirichlet(:u, Γ, (x, t) -> 0, Int[]) + @test_throws ErrorException("components not sorted: [2, 1]") Dirichlet(:u, Γ, (x, t) -> 0, Int[2, 1]) + @test_throws ErrorException("components not unique: [2, 2]") Dirichlet(:u, Γ, (x, t) -> 0, Int[2, 2]) + ## Scalar + dbc = Dirichlet(:s, Γ, (x, t) -> 0) + add!(ch, dbc) + @test dbc.components == [1] + dbc = Dirichlet(:s, Γ, (x, t) -> 0, [1]) + add!(ch, dbc) + @test dbc.components == [1] + dbc = Dirichlet(:s, Γ, (x, t) -> 0, [1, 2]) + @test_throws ErrorException("components [1, 2] not within range of field :s (1 dimension(s))") add!(ch, dbc) + dbc = Dirichlet(:p, Γ, (x, t) -> 0) + @test_throws ErrorException("could not find field :p in DofHandler (existing fields: [:s, :v])") add!(ch, dbc) + ## Vector + dbc = Dirichlet(:v, Γ, (x, t) -> 0) + add!(ch, dbc) + @test dbc.components == [1, 2] + dbc = Dirichlet(:v, Γ, (x, t) -> 0, [1]) + add!(ch, dbc) + @test dbc.components == [1] + dbc = Dirichlet(:v, Γ, (x, t) -> 0, [2, 3]) + @test_throws ErrorException("components [2, 3] not within range of field :v (2 dimension(s))") add!(ch, dbc) + + # PeriodicDirichlet + @test_throws ErrorException("components are empty: $(Int)[]") PeriodicDirichlet(:u, face_map, Int[]) + @test_throws ErrorException("components not sorted: [2, 1]") PeriodicDirichlet(:u, face_map, Int[2, 1]) + @test_throws ErrorException("components not unique: [2, 2]") PeriodicDirichlet(:u, face_map, Int[2, 2]) + ## Scalar + pdbc = PeriodicDirichlet(:s, face_map) + add!(ConstraintHandler(dh), pdbc) + @test pdbc.components == [1] + pdbc = PeriodicDirichlet(:s, face_map, (x,t) -> 0) + add!(ConstraintHandler(dh), pdbc) + @test pdbc.components == [1] + pdbc = PeriodicDirichlet(:s, face_map, [1]) + add!(ConstraintHandler(dh), pdbc) + @test pdbc.components == [1] + pdbc = PeriodicDirichlet(:s, face_map, [1, 2]) + @test_throws ErrorException("components [1, 2] not within range of field :s (1 dimension(s))") add!(ConstraintHandler(dh), pdbc) + pdbc = PeriodicDirichlet(:p, face_map) + @test_throws ErrorException("could not find field :p in DofHandler (existing fields: [:s, :v])") add!(ConstraintHandler(dh), pdbc) + ## Vector + pdbc = PeriodicDirichlet(:v, face_map) + add!(ConstraintHandler(dh), pdbc) + @test pdbc.components == [1, 2] + pdbc = PeriodicDirichlet(:v, face_map, (x, t) -> 0*x) + add!(ConstraintHandler(dh), pdbc) + @test pdbc.components == [1, 2] + pdbc = PeriodicDirichlet(:v, face_map, rand(2,2)) + add!(ConstraintHandler(dh), pdbc) + @test pdbc.components == [1, 2] + pdbc = PeriodicDirichlet(:v, face_map, rand(1,1)) + @test_throws ErrorException("size of rotation matrix does not match the number of components") add!(ConstraintHandler(dh), pdbc) + pdbc = PeriodicDirichlet(:v, face_map, (x, t) -> 0, [1]) + add!(ConstraintHandler(dh), pdbc) + @test pdbc.components == [1] + pdbc = PeriodicDirichlet(:v, face_map, (x, t) -> 0, [2, 3]) + @test_throws ErrorException("components [2, 3] not within range of field :v (2 dimension(s))") add!(ConstraintHandler(dh), pdbc) + pdbc = PeriodicDirichlet(:v, face_map, rand(2, 2), [2, 3]) + @test_throws ErrorException("components [2, 3] not within range of field :v (2 dimension(s))") add!(ConstraintHandler(dh), pdbc) +end + @testset "node bc" begin grid = generate_grid(Triangle, (1, 1)) addnodeset!(grid, "nodeset", x-> x[2] == -1 || x[1] == -1)