diff --git a/CHANGELOG.md b/CHANGELOG.md index 95a03f9732..ad2053acfe 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/download_resources.jl b/docs/download_resources.jl index 48b3a3600b..2c48d28848 100644 --- a/docs/download_resources.jl +++ b/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) diff --git a/docs/src/literate/computational_homogenization.jl b/docs/src/literate/computational_homogenization.jl index baa486f09c..33eabb607e 100644 --- a/docs/src/literate/computational_homogenization.jl +++ b/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") diff --git a/docs/src/literate/transient_heat.gif b/docs/src/literate/transient_heat.gif deleted file mode 100644 index 718c0d6ee0..0000000000 Binary files a/docs/src/literate/transient_heat.gif and /dev/null differ diff --git a/docs/src/literate/transient_heat_equation.jl b/docs/src/literate/transient_heat_equation.jl index c716b0af64..a5c354ef9c 100644 --- a/docs/src/literate/transient_heat_equation.jl +++ b/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, :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) @@ -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); diff --git a/docs/src/manual/boundary_conditions.md b/docs/src/manual/boundary_conditions.md index 0bcf5ec151..0a4de7b206 100644 --- a/docs/src/manual/boundary_conditions.md +++ b/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, :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{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}``. diff --git a/docs/src/reference/boundary_conditions.md b/docs/src/reference/boundary_conditions.md index 8f0f8c0a66..a87252f669 100644 --- a/docs/src/reference/boundary_conditions.md +++ b/docs/src/reference/boundary_conditions.md @@ -24,3 +24,9 @@ get_rhs_data apply_rhs! Ferrite.RHSData ``` + +# Initial conditions + +```@docs +apply_analytical! +``` diff --git a/src/Dofs/MixedDofHandler.jl b/src/Dofs/MixedDofHandler.jl index 85b91f599e..54ac744c70 100644 --- a/src/Dofs/MixedDofHandler.jl +++ b/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 diff --git a/src/Dofs/apply_analytical.jl b/src/Dofs/apply_analytical.jl new file mode 100644 index 0000000000..64099fd83d --- /dev/null +++ b/src/Dofs/apply_analytical.jl @@ -0,0 +1,98 @@ + +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, fieldname::Symbol, + f::Function, 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, fieldname::Symbol, f::Function, + 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, fieldname::Symbol, f::Function, + 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{dim,RefShape}, ip_geo::Interpolation, f::Function, cellset) where {dim, RefShape} + + coords = getcoordinates(dh.grid, first(cellset)) + ref_points = reference_coordinates(ip_fun) + dummy_weights = zeros(length(ref_points)) + qr = QuadratureRule{dim, RefShape}(dummy_weights, ref_points) + cv = CellScalarValues(qr, ip_fun, ip_geo) + 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 + getcoordinates!(coords, dh.grid, 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, cv, f) + end + return a +end + +function _apply_analytical!(a::Vector, dofs::Vector{Int}, coords::Vector{<:Vec}, field_dim, cv::CellScalarValues, f) + for i_dof in 1:getnquadpoints(cv) + x_dof = spatial_coordinate(cv, i_dof, coords) + for (idim, icval) in enumerate(f(x_dof)) + a[dofs[field_dim*(i_dof-1)+idim]] = icval + end + end + return a +end diff --git a/src/Ferrite.jl b/src/Ferrite.jl index 76d023abeb..63f26d8f46 100644 --- a/src/Ferrite.jl +++ b/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") diff --git a/src/Quadrature/quadrature.jl b/src/Quadrature/quadrature.jl index 71e56e4d53..a2c0de7375 100644 --- a/src/Quadrature/quadrature.jl +++ b/src/Quadrature/quadrature.jl @@ -40,6 +40,11 @@ struct QuadratureRule{dim,shape,T} points::Vector{Vec{dim,T}} end +function QuadratureRule{dim, shape}(weights::AbstractVector{Tw}, points::AbstractVector{Vec{dim,Tp}}) where {dim, shape, Tw, Tp} + T = promote_type(Tw, Tp) + QuadratureRule{dim,shape,T}(weights, points) +end + Base.copy(qr::QuadratureRule) = qr # TODO: Is it ever useful to get an actual copy? """ diff --git a/src/exports.jl b/src/exports.jl index 2f7b2cac70..9f3e8a3b75 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -111,6 +111,7 @@ export FieldHandler, Field, reshape_to_nodes, + apply_analytical!, # Constraints ConstraintHandler, diff --git a/test/runtests.jl b/test/runtests.jl index c9a2e3eaf1..edfb54ccf2 100644 --- a/test/runtests.jl +++ b/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 diff --git a/test/test_apply_analytical.jl b/test/test_apply_analytical.jl new file mode 100644 index 0000000000..b7136ca5a7 --- /dev/null +++ b/test/test_apply_analytical.jl @@ -0,0 +1,163 @@ +@testset "apply_analytical!" begin + + # Convenience helper functions + typebase(T::DataType) = T.name.wrapper + typebase(v) = typebase(typeof(v)) + change_ip_order(ip::Interpolation, ::Nothing) = ip + function change_ip_order(ip::Interpolation, order::Int) + B = typebase(ip) + Dim = Ferrite.getdim(ip) + RefShape = Ferrite.getrefshape(ip) + return B{Dim,RefShape,order}() + end + getcellorder(CT) = Ferrite.getorder(Ferrite.default_interpolation(CT)) + getcelltypedim(::Type{<:Cell{dim}}) where dim = dim + + # Functions to create dof handlers for testing + function testdh(CT, ip_order_u, ip_order_p) + dim = getcelltypedim(CT) + local grid + try + grid = generate_grid(CT, ntuple(_->3, dim)) + catch e + isa(e, MethodError) && e.f==generate_grid && return nothing + rethrow(e) + end + + dh = DofHandler(grid) + default_ip = Ferrite.default_interpolation(CT) + try + push!(dh, :u, dim, change_ip_order(default_ip, ip_order_u)) + push!(dh, :p, 1, change_ip_order(default_ip, ip_order_p)) + catch e + isa(e, MethodError) && e.f == Ferrite.reference_coordinates && return nothing + rethrow(e) + end + close!(dh) + return dh + end + + function testmdh(dim, ip_order_u, ip_order_p) + if dim == 1 + nodes = Node.([Vec{1}((x,)) for x in 0.0:1.0:3.0]) + cell1 = Cell{1,2,2}((1,2)) + cell2 = Cell{1,3,2}((2,3,4)) + grid = Grid([cell1, cell2], nodes; cellsets=Dict("A"=>Set(1:1), "B"=>Set(2:2))) + elseif dim == 2 + nodes = Node.([Vec{2}((x,y)) for y in 0.0:2 for x in 0.0:1]) + cell1 = Cell{2,3,3}((1,2,3)) + cell2 = Cell{2,3,3}((2,4,3)) + cell3 = Cell{2,4,4}((3,4,6,5)) + grid = Grid([cell1,cell2,cell3], nodes; cellsets=Dict("A"=>Set(1:2), "B"=>Set(3:3))) + else + error("Only dim=1 & 2 supported") + end + mdh = MixedDofHandler(grid) + default_ip_A = Ferrite.default_interpolation(getcelltype(grid, first(getcellset(grid,"A")))) + default_ip_B = Ferrite.default_interpolation(getcelltype(grid, first(getcellset(grid,"B")))) + ufield_A = Field(:u, change_ip_order(default_ip_A, ip_order_u), dim) + pfield_A = Field(:p, change_ip_order(default_ip_A, ip_order_p), 1) + ufield_B = Field(:u, change_ip_order(default_ip_B, ip_order_u), dim) + push!(mdh, FieldHandler([ufield_A, pfield_A], getcellset(grid,"A"))) + push!(mdh, FieldHandler([ufield_B,], getcellset(grid, "B"))) + close!(mdh) + return mdh + end + + # The following can be removed after #457 is merged if that will include the MixedDofHandler + function _global_dof_range(dh::MixedDofHandler, field_name::Symbol) + dofs = Set{Int}() + for fh in dh.fieldhandlers + if field_name ∈ Ferrite.getfieldnames(fh) + _global_dof_range!(dofs, dh, fh, field_name, fh.cellset) + end + end + return sort!(collect(Int, dofs)) + end + + function _global_dof_range(dh::DofHandler, field_name::Symbol) + dofs = Set{Int}() + _global_dof_range!(dofs, dh, dh, field_name, 1:getncells(dh.grid)) + return sort!(collect(Int, dofs)) + end + + function _global_dof_range!(dofs, dh, dh_fh, field_name, cellset) + eldofs = celldofs(dh, first(cellset)) + field_range = dof_range(dh_fh, field_name) + for i in cellset + celldofs!(eldofs, dh, i) + for j in field_range + @inbounds d = eldofs[j] + d in dofs || push!(dofs, d) + end + end + end + + @testset "DofHandler" begin + for (CT,name) in Ferrite.celltypes + for ip_order_u in 1:2 + for ip_order_p in 1:2 + dh = testdh(CT, ip_order_u, ip_order_p) + isnothing(dh) && continue # generate_grid not supported for this CT, or reference_coordinates not defined + dim = Ferrite.getdim(dh.grid) + num_udofs = length(_global_dof_range(dh, :u)) + num_pdofs = length(_global_dof_range(dh, :p)) + + # Test average value + a = zeros(ndofs(dh)) + f(x) = ones(Vec{dim}) + apply_analytical!(a, dh, :u, f) + @test sum(a)/length(a) ≈ num_udofs/(num_udofs+num_pdofs) + + # If not super/subparametric, compare with ConstraintHandler and node set + if ip_order_u==ip_order_p==getcellorder(CT) + fill!(a, 0) + a_ch = copy(a) + fp(x) = norm(x)^2 + fu(x::Vec{1}) = Vec{1}((sin(x[1]),)) + fu(x::Vec{2}) = Vec{2}((sin(x[1]), cos(x[2]))) + fu(x::Vec{3}) = Vec{3}((sin(x[1]), cos(x[2]), x[3]*(x[1]+x[2]))) + ch = ConstraintHandler(dh) + add!(ch, Dirichlet(:p, Set(1:getnnodes(dh.grid)), (x,t)->fp(x))) + add!(ch, Dirichlet(:u, Set(1:getnnodes(dh.grid)), (x,t)->fu(x))) + close!(ch); update!(ch, 0.0) + apply!(a_ch, ch) + + apply_analytical!(a, dh, :u, fu) + apply_analytical!(a, dh, :p, fp) + + @test a ≈ a_ch + end + end + end + end + end + + @testset "MixedDofHandler" begin + for dim in 1:2 + for ip_order_u in 1:2 + for ip_order_p in 1:2 + dh = testmdh(dim, ip_order_u, ip_order_p) + num_udofs = length(_global_dof_range(dh, :u)) + num_pdofs = length(_global_dof_range(dh, :p)) + + # Test average value + a = zeros(ndofs(dh)) + f(x) = ones(Vec{dim}) + apply_analytical!(a, dh, :u, f) + @test sum(a)/length(a) ≈ num_udofs/(num_udofs+num_pdofs) + end + end + end + end + + @testset "Exceptions" begin + dh = testdh(Quadrilateral, 1, 1) + @test_throws ErrorException apply_analytical!(zeros(ndofs(dh)), dh, :v, x->0.0) # Missing field + @test_throws ErrorException apply_analytical!(zeros(ndofs(dh)), dh, :u, x->0.0) # Should be f(x)::Vec{2} + + mdh = testmdh(2, 1, 1) + @test_throws ErrorException apply_analytical!(zeros(ndofs(mdh)), mdh, :v, x->0.0) # Missing field + @test_throws ErrorException apply_analytical!(zeros(ndofs(mdh)), mdh, :u, x->0.0) # Should be f(x)::Vec{2} + end +end