Add apply_analytical! function

CHANGELOG.md

@@ -12,6 +12,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
    support a new DoF renumbering order `DofOrder.Ext{Metis}()` that can be passed to
    `renumber!` to renumber DoFs using the Metis.jl library. ([#393][github-393],
    [#549][github-549])
+ - New function `apply_analytical!` for setting the values of the degrees of freedom for a 
+   specific field according to a spatial function `f(x)`. ([#532][github-532])
+
 
 ## [0.3.10] - 2022-12-11
 ### Added
@@ -210,6 +213,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 [github-528]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/528
 [github-529]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/529
 [github-530]: https://github.com/Ferrite-FEM/Ferrite.jl/issues/530
+[github-532]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/532
 [github-533]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/533
 [github-534]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/534
 [github-535]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/535

docs/download_resources.jl

@@ -7,6 +7,8 @@ mkpath(directory)
 for (file, url) in [
         "periodic-rve.msh" => "https://raw.githubusercontent.com/Ferrite-FEM/Ferrite.jl/gh-pages/assets/periodic-rve.msh",
         "periodic-rve-coarse.msh" => "https://raw.githubusercontent.com/Ferrite-FEM/Ferrite.jl/gh-pages/assets/periodic-rve-coarse.msh",
+        "transient_heat.gif" => "https://raw.githubusercontent.com/Ferrite-FEM/Ferrite.jl/gh-pages/assets/transient_heat.gif",
+        "transient_heat_colorbar.svg" => "https://raw.githubusercontent.com/Ferrite-FEM/Ferrite.jl/gh-pages/assets/transient_heat_colorbar.svg",
     ]
     afile = joinpath(directory, file)
     if !isfile(afile)

docs/src/literate/computational_homogenization.jl

@@ -516,16 +516,12 @@ round.(ev; digits=-8)
 # Finally, we export the solution and the stress field to a VTK file. For the export we
 # also compute the macroscopic part of the displacement.
 
-chM = ConstraintHandler(dh)
-add!(chM, Dirichlet(:u, Set(1:getnnodes(grid)), (x, t) -> εᴹ[Int(t)] ⋅ x, [1, 2]))
-close!(chM)
 uM = zeros(ndofs(dh))
 
 vtk_grid("homogenization", dh) do vtk
     for i in 1:3
         ## Compute macroscopic solution
-        update!(chM, i)
-        apply!(uM, chM)
+        apply_analytical!(uM, dh, x -> εᴹ[i] ⋅ x, :u)
         ## Dirichlet
         vtk_point_data(vtk, dh, uM + u.dirichlet[i], "_dirichlet_$i")
         vtk_point_data(vtk, projector, σ.dirichlet[i], "σvM_dirichlet_$i")

docs/src/literate/transient_heat_equation.jl

@@ -1,6 +1,7 @@
 # # Time Dependent Problems
 #
 # ![](transient_heat.gif)
+# ![](transient_heat_colorbar.svg)
 #
 # *Figure 1*: Visualization of the temperature time evolution on a unit
 # square where the prescribed temperature on the upper and lower parts
@@ -20,8 +21,8 @@
 # ```
 #
 # where $u$ is the unknown temperature field, $k$ the heat conductivity,
-# $f$ the heat source and $\Omega$ the domain. For simplicity we set $f = 1$
-# and $k = 1$. We define homogeneous Dirichlet boundary conditions along the left and right edge of the domain.
+# $f$ the heat source and $\Omega$ the domain. For simplicity, we hard code $f = 0.1$
+# and $k = 10^{-3}$. We define homogeneous Dirichlet boundary conditions along the left and right edge of the domain.
 # ```math
 # u(x,t) = 0 \quad x \in \partial \Omega_1,
 # ```
@@ -34,12 +35,12 @@
 # ```
 # The semidiscrete weak form is given by
 # ```math
-# \int_{\Omega}v \frac{\partial u}{\partial t} \ \mathrm{d}\Omega + \int_{\Omega} \nabla v \cdot \nabla u \ \mathrm{d}\Omega = \int_{\Omega} v \ \mathrm{d}\Omega,
+# \int_{\Omega}v \frac{\partial u}{\partial t} \ \mathrm{d}\Omega + \int_{\Omega} k \nabla v \cdot \nabla u \ \mathrm{d}\Omega = \int_{\Omega} f v \ \mathrm{d}\Omega,
 # ```
 # where $v$ is a suitable test function. Now, we still need to discretize the time derivative. An implicit Euler scheme is applied,
 # which yields:
 # ```math
-# \int_{\Omega} v\, u_{n+1}\ \mathrm{d}\Omega + \Delta t\int_{\Omega} \nabla v \cdot \nabla u_{n+1} \ \mathrm{d}\Omega = \Delta t\int_{\Omega} v \ \mathrm{d}\Omega + \int_{\Omega} v \, u_{n} \ \mathrm{d}\Omega.
+# \int_{\Omega} v\, u_{n+1}\ \mathrm{d}\Omega + \Delta t\int_{\Omega} k \nabla v \cdot \nabla u_{n+1} \ \mathrm{d}\Omega = \Delta t\int_{\Omega} f v \ \mathrm{d}\Omega + \int_{\Omega} v \, u_{n} \ \mathrm{d}\Omega.
 # ```
 # If we assemble the discrete operators, we get the following algebraic system:
 # ```math
@@ -91,15 +92,17 @@ f = zeros(ndofs(dh));
 max_temp = 100
 Δt = 1
 T = 200
+t_rise = 100
 ch = ConstraintHandler(dh);
 
 # Here, we define the boundary condition related to $\partial \Omega_1$.
-∂Ω₁ = union(getfaceset.((grid, ), ["left", "right"])...)
+∂Ω₁ = union(getfaceset.((grid,), ["left", "right"])...)
 dbc = Dirichlet(:u, ∂Ω₁, (x, t) -> 0)
 add!(ch, dbc);
-# While the next code block corresponds to the linearly increasing temperature description on $\partial \Omega_2$.
-∂Ω₂ = union(getfaceset.((grid, ), ["top", "bottom"])...)
-dbc = Dirichlet(:u, ∂Ω₂, (x, t) -> t*(max_temp/T))
+# While the next code block corresponds to the linearly increasing temperature description on $\partial \Omega_2$
+# until `t=t_rise`, and then keep constant
+∂Ω₂ = union(getfaceset.((grid,), ["top", "bottom"])...)
+dbc = Dirichlet(:u, ∂Ω₂, (x, t) -> max_temp * clamp(t / t_rise, 0, 1))
 add!(ch, dbc)
 close!(ch)
 update!(ch, 0.0);
@@ -125,12 +128,12 @@ function doassemble_K!(K::SparseMatrixCSC, f::Vector, cellvalues::CellScalarValu
             dΩ = getdetJdV(cellvalues, q_point)
 
             for i in 1:n_basefuncs
-                v  = shape_value(cellvalues, q_point, i)
+                v = shape_value(cellvalues, q_point, i)
                 ∇v = shape_gradient(cellvalues, q_point, i)
-                fe[i] += v * dΩ
+                fe[i] += 0.1 * v * dΩ
                 for j in 1:n_basefuncs
                     ∇u = shape_gradient(cellvalues, q_point, j)
-                    Ke[i, j] += (∇v ⋅ ∇u) * dΩ
+                    Ke[i, j] += 1e-3 * (∇v ⋅ ∇u) * dΩ
                 end
             end
         end
@@ -158,7 +161,7 @@ function doassemble_M!(M::SparseMatrixCSC, cellvalues::CellScalarValues{dim}, dh
             dΩ = getdetJdV(cellvalues, q_point)
 
             for i in 1:n_basefuncs
-                v  = shape_value(cellvalues, q_point, i)
+                v = shape_value(cellvalues, q_point, i)
                 for j in 1:n_basefuncs
                     u = shape_value(cellvalues, q_point, j)
                     Me[i, j] += (v * u) * dΩ
@@ -180,16 +183,24 @@ A = (Δt .* K) + M;
 # by `get_rhs_data`. The function returns a `RHSData` struct, which contains all needed informations to apply
 # the boundary conditions solely on the right-hand-side vector of the problem.
 rhsdata = get_rhs_data(ch, A);
-# We set the initial time step, denoted by uₙ,  to $\mathbf{0}$.
+# We set the values at initial time step, denoted by uₙ, to a bubble-shape described by 
+# $(x_1^2-1)(x_2^2-1)$, such that it is zero at the boundaries and the maximum temperature in the center.
 uₙ = zeros(length(f));
+apply_analytical!(uₙ, dh, x -> (x[1]^2 - 1) * (x[2]^2 - 1) * max_temp, :u);
 # Here, we apply **once** the boundary conditions to the system matrix `A`.
 apply!(A, ch);
 
 # To store the solution, we initialize a `paraview_collection` (.pvd) file.
 pvd = paraview_collection("transient-heat.pvd");
+t = 0
+vtk_grid("transient-heat-$t", dh) do vtk
+    vtk_point_data(vtk, dh, uₙ)
+    vtk_save(vtk)
+    pvd[t] = vtk
+end
 
 # At this point everything is set up and we can finally approach the time loop.
-for t in 0:Δt:T
+for t in Δt:Δt:T
     #First of all, we need to update the Dirichlet boundary condition values.
     update!(ch, t)
 
@@ -199,15 +210,15 @@ for t in 0:Δt:T
     apply_rhs!(rhsdata, b, ch)
 
     #Finally, we can solve the time step and save the solution afterwards.
-    u = A \ b;
+    u = A \ b
 
     vtk_grid("transient-heat-$t", dh) do vtk
         vtk_point_data(vtk, dh, u)
         vtk_save(vtk)
         pvd[t] = vtk
     end
-   #At the end of the time loop, we set the previous solution to the current one and go to the next time step.
-   uₙ .= u
+    #At the end of the time loop, we set the previous solution to the current one and go to the next time step.
+    uₙ .= u
 end
 # In order to use the .pvd file we need to store it to the disk, which is done by:
 vtk_save(pvd);

docs/src/manual/boundary_conditions.md

@@ -260,3 +260,34 @@ pdbc = PeriodicDirichlet(
         (x, t) -> ū + ∇ū  ⋅ (x - x̄)
     )
     ```
+
+# Initial Conditions
+
+When solving time-dependent problems, initial conditions, different from zero, may be required. 
+For finite element formulations of ODE-type, 
+i.e. ``\boldsymbol{u}'(t) = \boldsymbol{f}(\boldsymbol{u}(t),t)``, 
+where ``\boldsymbol{u}(t)`` are the degrees of freedom,
+initial conditions can be specified by the [`apply_analytical!`](@ref) function.
+For example, specify the initial pressure as a function of the y-coordinate
+```julia
+ρ = 1000; g = 9.81    # density [kg/m³] and gravity [N/kg]
+grid = generate_grid(Quadrilateral, (10,10))
+dh = DofHandler(grid); push!(dh, :u, 2); push!(dh, :p, 1); close!(dh)
+u = zeros(ndofs(dh))
+apply_analytical!(u, dh, x -> ρ * g * x[2], :p)
+```
+
+See also [Time Dependent Problems](@ref) for one example. 
+
+*Note about solving DAE:* 
+A Differential Algebraic Equations (DAE) is an equation of the form
+``\boldsymbol{r}(\boldsymbol{u}(t),\boldsymbol{u}'(t),t)=\boldsymbol{0}``,
+which usually cannot be expressed as a true ODE. They occur often,
+but not always, in forms where some time derivatives are missing
+```math
+u_1'(t) = f(\boldsymbol{u}(t),t)
+0 = g(\boldsymbol{u}(t),t)`
+```
+In for such equations, it is usually necessary to specify initial conditions 
+for both ``\boldsymbol{u}(0)`` and ``\boldsymbol{u}'(0)``, and these must be consistent,
+i.e. ``\boldsymbol{r}(\boldsymbol{u}(0),\boldsymbol{u}'(0),0)=\boldsymbol{0}``.

docs/src/reference/boundary_conditions.md

@@ -24,3 +24,9 @@ get_rhs_data
 apply_rhs!
 Ferrite.RHSData
 ```
+
+# Initial conditions
+
+```@docs
+apply_analytical!
+```

src/Dofs/MixedDofHandler.jl

@@ -400,7 +400,9 @@ end
 
 find_field(dh::MixedDofHandler, field_name::Symbol) = find_field(first(dh.fieldhandlers), field_name)
 field_offset(dh::MixedDofHandler, field_name::Symbol) = field_offset(first(dh.fieldhandlers), field_name)
+getfieldinterpolation(fh::FieldHandler, field_idx::Int) = fh.fields[field_idx].interpolation
 getfieldinterpolation(dh::MixedDofHandler, field_idx::Int) = dh.fieldhandlers[1].fields[field_idx].interpolation
+getfielddim(fh::FieldHandler, field_idx::Int) = fh.fields[field_idx].dim
 getfielddim(dh::MixedDofHandler, field_idx::Int) = dh.fieldhandlers[1].fields[field_idx].dim
 
 function reshape_to_nodes(dh::MixedDofHandler, u::Vector{T}, fieldname::Symbol) where T

src/Dofs/apply_analytical.jl

@@ -0,0 +1,102 @@
+
+function _default_interpolations(dh::MixedDofHandler)
+    fhs = dh.fieldhandlers
+    getcelltype(i) = typeof(getcells(dh.grid, first(fhs[i].cellset)))
+    ntuple(i -> default_interpolation(getcelltype(i)), length(fhs))
+end
+
+function _default_interpolation(dh::DofHandler)
+    return default_interpolation(typeof(getcells(dh.grid, 1)))
+end
+
+"""
+    apply_analytical!(
+        a::AbstractVector, dh::AbstractDofHandler, f::Function, 
+        fieldname::Symbol, cellset=1:getncells(dh.grid))
+
+Apply a solution `f(x)` by modifying the values in the degree of freedom vector `a`
+pertaining to the field `fieldname` for all cells in `cellset`.
+The function `f(x)` are given the spatial coordinate 
+of the degree of freedom. For scalar fields, `f(x)::Number`, 
+and for vector fields with dimension `dim`, `f(x)::Vec{dim}`.
+
+This function can be used to apply initial conditions for time dependent problems.
+
+!!! note
+    This function only works for standard nodal finite element interpolations 
+    when the function value at the (algebraic) node is equal to the corresponding 
+    degree of freedom value.
+    This holds for e.g. Lagrange and Serendipity interpolations, including 
+    sub- and superparametric elements. 
+
+"""
+function apply_analytical!(
+    a::AbstractVector, dh::DofHandler, f::Function, fieldname::Symbol,
+    cellset=1:getncells(dh.grid)
+    )
+    fieldname ∉ getfieldnames(dh) && error("The fieldname $fieldname was not found in the dof handler")
+    ip_geo = _default_interpolation(dh)
+    field_idx = find_field(dh, fieldname)
+    ip_fun = getfieldinterpolation(dh, field_idx)
+    celldofinds = dof_range(dh, fieldname)
+    field_dim = getfielddim(dh, field_idx)
+    _apply_analytical!(a, dh, celldofinds, field_dim, ip_fun, ip_geo, f, cellset)
+end
+
+function apply_analytical!(
+    a::AbstractVector, dh::MixedDofHandler, f::Function, fieldname::Symbol,
+    cellset=1:getncells(dh.grid),
+    )
+    fieldname ∉ getfieldnames(dh) && error("The fieldname $fieldname was not found in the dof handler")
+    ip_geos = _default_interpolations(dh)
+
+    for (fh, ip_geo) in zip(dh.fieldhandlers, ip_geos)
+        fieldname ∈ getfieldnames(fh) || continue
+        field_idx = find_field(fh, fieldname)
+        ip_fun = getfieldinterpolation(fh, field_idx)
+        field_dim = getfielddim(fh, field_idx)
+        celldofinds = dof_range(fh, fieldname)
+        _apply_analytical!(a, dh, celldofinds, field_dim, ip_fun, ip_geo, f, intersect(fh.cellset, cellset))
+    end
+    return a
+end
+
+function _apply_analytical!(
+    a::Vector, dh::AbstractDofHandler, celldofinds, field_dim,
+    ip_fun::Interpolation, ip_geo::Interpolation, f::Function, cellset)
+
+    coords = getcoordinates(dh.grid, first(cellset))
+    c_dofs = celldofs(dh, first(cellset))
+    f_dofs = zeros(Int, length(celldofinds))
+
+    # Check f before looping
+    length(f(first(coords))) == field_dim || error("length(f(x)) must be equal to dimension of the field ($field_dim)")
+
+    for cellnr in cellset
+        cellcoords!(coords, dh, cellnr)
+        celldofs!(c_dofs, dh, cellnr)
+        for (i, celldofind) in enumerate(celldofinds)
+            f_dofs[i] = c_dofs[celldofind]
+        end
+        _apply_analytical!(a, f_dofs, coords, field_dim, ip_fun, ip_geo, f)
+    end
+    return a
+end
+
+function _apply_analytical!(a::Vector, dofs::Vector{Int}, coords::Vector{XT}, field_dim, ip_fun, ip_geo, f) where {XT<:Vec}
+
+    getnbasefunctions(ip_geo) == length(coords) || error("coords=$coords not compatible with ip_ge=$ip_geo")
+    ref_coords = reference_coordinates(ip_fun)
+    length(ref_coords) * field_dim == length(dofs) || error("$ip_fun must have length(dofs)=$(length(dofs)) reference coords")
+
+    for (i_dof, refpoint) in enumerate(ref_coords)
+        x_dof = zero(XT)
+        for (i, x_node) in enumerate(coords)
+            x_dof += value(ip_geo, i, refpoint) * x_node
+        end
+        for (idim, icval) in enumerate(f(x_dof))
+            a[dofs[field_dim*(i_dof-1)+idim]] = icval
+        end
+    end
+    return a
+end

src/Ferrite.jl

@@ -60,6 +60,7 @@ include("Grid/coloring.jl")
 include("Dofs/DofHandler.jl")
 include("Dofs/MixedDofHandler.jl")
 include("Dofs/ConstraintHandler.jl")
+include("Dofs/apply_analytical.jl")
 include("Dofs/DofRenumbering.jl")
 
 include("iterators.jl")

src/exports.jl

@@ -111,6 +111,7 @@ export
     FieldHandler,
     Field,
     reshape_to_nodes,
+    apply_analytical!,
 
 # Constraints
     ConstraintHandler,

test/runtests.jl

@@ -27,5 +27,6 @@ include("test_l2_projection.jl")
 include("test_pointevaluation.jl")
 # include("test_notebooks.jl")
 include("test_apply_rhs.jl")
+include("test_apply_analytical.jl")
 include("test_examples.jl")
 @test all(x -> isdefined(Ferrite, x), names(Ferrite))  # Test that all exported symbols are defined