diff --git a/docs/src/literate/heat_equation.jl b/docs/src/literate/heat_equation.jl index 118cd1074c..0095901cc6 100644 --- a/docs/src/literate/heat_equation.jl +++ b/docs/src/literate/heat_equation.jl @@ -93,19 +93,22 @@ ch = ConstraintHandler(dh); # Now we are set up to define our constraint. We specify which field # the condition is for, and our combined face set `∂Ω`. The last -# argument is a function which takes the spatial coordinate $\textbf{x}$ and -# the current time $t$ and returns the prescribed value. In this case -# it is trivial -- no matter what $\textbf{x}$ and $t$ we return $0$. When we have +# argument is a function of the form $f(\textbf{x})$ or $f(\textbf{x}, t)$, +# where $\textbf{x}$ is the spatial coordinate and +# $t$ the current time, and returns the prescribed value. Since the boundary condition in +# this case do not depend on time we define our function as $f(\textbf{x}) = 0$, i.e. +# no matter what $\textbf{x}$ we return $0$. When we have # specified our constraint we `add!` it to `ch`. dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0) add!(ch, dbc); -# We also need to `close!` and `update!` our boundary conditions. When we call `close!` +# Finally we also need to `close!` our constraint handler. When we call `close!` # the dofs corresponding to our constraints are calculated and stored -# in our `ch` object. Since the boundary conditions are, in this case, -# independent of time we can `update!` them directly with e.g. $t = 0$. +# in our `ch` object. close!(ch) -update!(ch, 0.0); + +# Note that if one or more of the constraints are time dependent we would use +# [`update!`](@ref) to recompute prescribed values in each new timestep. # ### Assembling the linear system # diff --git a/docs/src/manual/boundary_conditions.md b/docs/src/manual/boundary_conditions.md index 0a4de7b206..f4cd1a49d4 100644 --- a/docs/src/manual/boundary_conditions.md +++ b/docs/src/manual/boundary_conditions.md @@ -27,14 +27,14 @@ for computing the prescribed value. Example: dbc1 = Dirichlet( :u, # Name of the field getfaceset(grid, "left"), # Part of the boundary - (x, t) -> 1.0 * t, # Function mapping coordinate and time to a prescribed value + x -> 1.0, # Function mapping coordinate to a prescribed value ) ``` The field name is given as a symbol, just like when the field was added to the dof handler, the part of the boundary where this constraint is active is given as a face set, and the -function computing the prescribed value should accept two input arguments (coordinate `x` -and time `t`) and return the prescribed value. +function computing the prescribed value should be of the form `f(x)` or `f(x, t)` +(coordinate `x` and time `t`) and return the prescribed value(s). !!! note "Multiple sets" To apply a constraint on multiple face sets in the grid you can use `union` to join @@ -45,8 +45,8 @@ and time `t`) and return the prescribed value. creates a new face set containing all faces in the `"left"` and "`right`" face sets, which can be passed to the `Dirichlet` constructor. -By default the constraint is added to the first component of the given field. To add the -constraint to multiple components a fourth argument with the components should be passed to +By default the constraint is added to all components of the given field. To add the +constraint to selected components a fourth argument with the components should be passed to the constructor. Here is an example where a constraint is added to component 1 and 3 of a vector field `:u`: @@ -54,7 +54,7 @@ vector field `:u`: dbc2 = Dirichlet( :u, # Name of the field getfaceset(grid, "left"), # Part of the boundary - (x, t) -> [0.0, 0.0], # Function mapping coordinate and time to a prescribed value + x -> [0.0, 0.0], # Function mapping coordinate to prescribed values [1, 3], # Components ) ``` @@ -73,10 +73,9 @@ Finally, just like for the dof handler, we need to use [`close!`](@ref) to final constraint handler. Internally this will then compute the degrees-of-freedom that match the constraints we added. -Since the constraints can in general depend on time we also need to need to call -[`update!`](@ref) with the current time in order to compute the prescribed values. The -same constraint handler can then be used for all time steps by calling `update!` with the -proper time, e.g.: +If one or more of the constraints depend on time, i.e. they are specified as `f(x, t)`, the +prescribed values can be recomputed in each new time step by calling [`update!`](@ref) with +the proper time, e.g.: ```julia for t in 0.0:0.1:1.0 @@ -85,10 +84,6 @@ for t in 0.0:0.1:1.0 end ``` -!!! note - You *must* call `update!`, even if your constraints does not depend on time - (as `dbc2` above), e.g. `update!(ch, 0.0)`. - !!! note "Examples" Most examples make use of Dirichlet boundary conditions, for example [Heat Equation](@ref). diff --git a/docs/src/reference/boundary_conditions.md b/docs/src/reference/boundary_conditions.md index a87252f669..950f3fd2d7 100644 --- a/docs/src/reference/boundary_conditions.md +++ b/docs/src/reference/boundary_conditions.md @@ -16,6 +16,7 @@ collect_periodic_faces collect_periodic_faces! add! close! +update! apply! apply_zero! apply_local! diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index d4f941bc94..6b512bc9a0 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -3,7 +3,7 @@ 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` +the boundary. `f` is a function of the form `f(x)` or `f(x, t)` where `x` is the spatial coordinate and `t` is the current time, and returns the prescribed value. `components` specify the components of `u` that are prescribed by this condition. By default all components @@ -22,7 +22,7 @@ Dirichlet condition for the `:u` field, on the faceset called dbc = Dirichlet(:s, ∂Ω, (x, t) -> sin(t)) # Prescribe all components of vector field :v on ∂Ω to 0 -dbc = Dirichlet(:v, ∂Ω, (x, t) -> 0 * x) +dbc = Dirichlet(:v, ∂Ω, x -> 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]) @@ -32,7 +32,7 @@ dbc = Dirichlet(:v, ∂Ω, (x, t) -> [sin(t), cos(t)], [2, 3]) which applies the condition via [`apply!`](@ref) and/or [`apply_zero!`](@ref). """ struct Dirichlet # <: Constraint - f::Function # f(x,t) -> value + f::Function # f(x) or f(x,t) -> value(s) faces::Union{Set{Int},Set{FaceIndex},Set{EdgeIndex},Set{VertexIndex}} field_name::Symbol components::Vector{Int} # components of the field @@ -222,6 +222,12 @@ function close!(ch::ConstraintHandler) end ch.closed[] = true + + # Compute the prescribed values by calling update!: This should be cheap, and for the + # common case where constraints does not depend on time it is annoying and easy to + # forget to call this on the outside. + update!(ch) + return ch end @@ -392,13 +398,25 @@ function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::Set{Int}, interpo return ch end -# Updates the DBC's to the current time `time` +""" + update!(ch::ConstraintHandler, time::Real=0.0) + +Update time-dependent inhomogeneities for the new time. This calls `f(x)` or `f(x, t)` when +applicable, where `f` is the function(s) corresponding to the constraints in the handler, to +compute the inhomogeneities. + +Note that this is called implicitly in `close!(::ConstraintHandler)`. +""" function update!(ch::ConstraintHandler, time::Real=0.0) @assert ch.closed[] - for (i,dbc) in enumerate(ch.dbcs) + for (i, dbc) in pairs(ch.dbcs) + # If the BC function only accept one argument, i.e. f(x), we create a wrapper + # g(x, t) = f(x) that discards the second parameter so that _update! can always call + # the function with two arguments internally. + wrapper_f = hasmethod(dbc.f, Tuple{Any,Any}) ? dbc.f : (x, _) -> dbc.f(x) # Function barrier - _update!(ch.inhomogeneities, dbc.f, dbc.faces, dbc.field_name, dbc.local_face_dofs, dbc.local_face_dofs_offset, - dbc.components, ch.dh, ch.bcvalues[i], ch.dofmapping, ch.dofcoefficients, convert(Float64, time)) + _update!(ch.inhomogeneities, wrapper_f, dbc.faces, dbc.field_name, dbc.local_face_dofs, dbc.local_face_dofs_offset, + dbc.components, ch.dh, ch.bcvalues[i], ch.dofmapping, ch.dofcoefficients, time) end # Compute effective inhomogeneity for affine constraints with prescribed dofs in the # RHS. For example, in u2 = w3 * u3 + w4 * u4 + b2 we allow e.g. u3 to be prescribed by @@ -425,7 +443,7 @@ end # for faces function _update!(inhomogeneities::Vector{Float64}, f::Function, faces::Set{<:BoundaryIndex}, field::Symbol, local_face_dofs::Vector{Int}, local_face_dofs_offset::Vector{Int}, components::Vector{Int}, dh::AbstractDofHandler, facevalues::BCValues, - dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::T) where {T} + dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::Real) where {T} cc = CellCache(dh, UpdateFlags(; nodes=false, coords=true, dofs=true)) for (cellidx, faceidx) in faces @@ -463,7 +481,7 @@ end # for nodes function _update!(inhomogeneities::Vector{Float64}, f::Function, nodes::Set{Int}, field::Symbol, nodeidxs::Vector{Int}, globaldofs::Vector{Int}, components::Vector{Int}, dh::AbstractDofHandler, facevalues::BCValues, - dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::T) where T + dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::Real) where T counter = 1 for (idx, nodenumber) in enumerate(nodeidxs) x = dh.grid.nodes[nodenumber].x diff --git a/test/test_constraints.jl b/test/test_constraints.jl index 61dc63d2eb..e108f7b914 100644 --- a/test/test_constraints.jl +++ b/test/test_constraints.jl @@ -1007,8 +1007,8 @@ end # testset compare_by_dbc( dh, PeriodicDirichlet(:v, collect_periodic_faces(grid, "left", "right"), collect(1:D)), - Dirichlet(:v, getfaceset(grid, "left"), (x,t) -> [0., 0.], collect(1:D)), - Dirichlet(:v, getfaceset(grid, "right"), (x,t) -> [0., 0.], collect(1:D)), + Dirichlet(:v, getfaceset(grid, "left"), (x,t) -> fill(0., D), collect(1:D)), + Dirichlet(:v, getfaceset(grid, "right"), (x,t) -> fill(0., D), collect(1:D)), ) compare_by_dbc( dh, @@ -1032,8 +1032,8 @@ end # testset compare_by_dbc( dh, PeriodicDirichlet(:v, collect_periodic_faces(grid, "front", "back"), 1:D), - Dirichlet(:v, getfaceset(grid, "front"), (x,t) -> [0., 0.], 1:D), - Dirichlet(:v, getfaceset(grid, "back"), (x,t) -> [0., 0.], 1:D), + Dirichlet(:v, getfaceset(grid, "front"), (x,t) -> fill(0., D), 1:D), + Dirichlet(:v, getfaceset(grid, "back"), (x,t) -> fill(0., D), 1:D), ) compare_by_dbc( dh,