Add apply_analytical! function

docs/generate.jl

@@ -41,7 +41,7 @@ include("download_resources.jl")
                 Literate.notebook(input, GENERATEDDIR, preprocess = nbpre, execute = is_ci) # Don't execute locally
             end
         end
-    elseif any(endswith.(example, [".png", ".jpg", ".gif"]))
+    elseif any(endswith.(example, [".png", ".jpg", ".gif", ".svg"]))
         cp(joinpath(EXAMPLEDIR, example), joinpath(GENERATEDDIR, example); force=true)
     else
         @warn "ignoring $example"

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, :u, x->εᴹ[i] ⋅ x)
         ## 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"])...)
 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$.
+# While the next code block corresponds to the linearly increasing temperature description on $\partial \Omega_2$
+# untill `t=t_rise`, and then keep constant
 ∂Ω₂ = union(getfaceset.((grid, ), ["top", "bottom"])...)
-dbc = Dirichlet(:u, ∂Ω₂, (x, t) -> t*(max_temp/T))
+dbc = Dirichlet(:u, ∂Ω₂, (x, t) -> max_temp*clamp(t/t_rise, 0, 1))
 add!(ch, dbc)
 close!(ch)
 update!(ch, 0.0);
@@ -127,10 +130,10 @@ function doassemble_K!(K::SparseMatrixCSC, f::Vector, cellvalues::CellScalarValu
             for i in 1:n_basefuncs
                 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] += 1.e-3*(∇v ⋅ ∇u) * dΩ
                 end
             end
         end
@@ -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 half the max temperature in the center.
 uₙ = zeros(length(f));
+apply_analytical!(uₙ, dh, :u, x->(x[1]^2-1)*(x[2]^2-1)*max_temp);
 # 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)
 

docs/src/manual/boundary_conditions.md

@@ -260,3 +260,29 @@ 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{x}'(t) = \boldsymbol{f}(\boldsymbol{x}(t),t)``, 
+where ``\boldsymbol{x}(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)
+x = zeros(ndofs(dh))
+apply_analytical!(x, dh, :p, x->ρ*g*x[2])
+```
+
+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{x}(t),\boldsymbol{x}'(t),t)=\boldsymbol{0}``,
+which cannot be expressed as an ODE.
+In for such equations, it is usually necessary to specify initial conditions 
+for both ``\boldsymbol{x}(0)`` and ``\boldsymbol{x}'(0)``, and these must be consistent,
+i.e. ``\boldsymbol{r}(\boldsymbol{x}(0),\boldsymbol{x}'(0),0)=\boldsymbol{0}``.

docs/src/reference/boundary_conditions.md

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

src/Dofs/MixedDofHandler.jl

@@ -463,7 +463,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/Ferrite.jl

@@ -59,6 +59,7 @@ include("Dofs/MixedDofHandler.jl")
 include("Dofs/ConstraintHandler.jl")
 
 include("iterators.jl")
+include("initial_conditions.jl")
 
 # Assembly
 include("assembler.jl")

src/exports.jl

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