Implicitly call update!(ch) in close!(ch) and allow f(x). (#459)

This patch adds an implicit call to update! when closing the
ConstraintHandler. This step is easy to forget, in particular if your
problem doesn't depend on time at all. In addition, it is now possible
to specify constraint functions of the form f(x) directly.

Fixes #207, closes #213, closes #435.

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
 #

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,16 +45,16 @@ 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`:
 
 ```julia
 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).

docs/src/reference/boundary_conditions.md

@@ -16,6 +16,7 @@ collect_periodic_faces
 collect_periodic_faces!
 add!
 close!
+update!
 apply!
 apply_zero!
 apply_local!

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

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,