Increased test coverage for internal operations (#104)

src/utils.jl

@@ -512,17 +512,15 @@ function interpolate_gradient_field(dh::Ferrite.DofHandler{spatial_dim}, u::Abst
     dh_gradient = Ferrite.DofHandler(getgrid(dh))
     ip_gradient = get_gradient_interpolation(ip)
     field_dim = Ferrite.getfielddim(dh,field_name)
-    push!(dh_gradient, :gradient, field_dim*spatial_dim, ip_gradient) # field dim × spatial dim components
+    Ferrite.add!(dh_gradient, :gradient, field_dim*spatial_dim, ip_gradient) # field dim × spatial dim components
     for fieldname in copy_fields
         _field_idx = Ferrite.find_field(dh, fieldname)
         _ip = Ferrite.getfieldinterpolation(dh, _field_idx)
         _field_dim = Ferrite.getfielddim(dh,fieldname)
-        push!(dh_gradient, fieldname, _field_dim, _ip)
+        Ferrite.add!(dh_gradient, fieldname, _field_dim, _ip)
     end
     Ferrite.close!(dh_gradient)
 
-    num_base_funs = Ferrite.getnbasefunctions(ip_gradient)
-
     # FIXME this does not work for mixed grids
     ip_geom = Ferrite.default_interpolation(typeof(Ferrite.getcells(getgrid(dh), 1)))
     ref_coords_gradient = Ferrite.reference_coordinates(ip_gradient)
@@ -535,22 +533,22 @@ function interpolate_gradient_field(dh::Ferrite.DofHandler{spatial_dim}, u::Abst
 
     # Allocate storage for the fluxes to store
     u_gradient = zeros(Ferrite.ndofs(dh_gradient))
-    # In general uᵉ_gradient is an order 3 tensor [field_dim, spatial_dim, num_base_funs]
+    # In general uᵉ_gradient is an order 3 tensor [field_dim, spatial_dim, nqp]
     uᵉ_gradient = zeros(length(cell_dofs_gradient[Ferrite.dof_range(dh_gradient, :gradient)]))
-    uᵉ_gradient_view = reshape(uᵉ_gradient, (spatial_dim, field_dim, num_base_funs))
-    uᵉ = zeros(field_dim*Ferrite.getnbasefunctions(ip))
+    uᵉshape = (spatial_dim, field_dim, Ferrite.getnquadpoints(cv))
+    uᵉ_gradient_view = reshape(uᵉ_gradient, uᵉshape)
 
     for (cell_num, cell) in enumerate(Ferrite.CellIterator(dh))
         # Get element dofs on parent field
         Ferrite.celldofs!(cell_dofs, dh, cell_num)
-        uᵉ .= u[cell_dofs[Ferrite.dof_range(dh, field_name)]]
+        uᵉ = @views u[cell_dofs[Ferrite.dof_range(dh, field_name)]]
 
         # And initialize cellvalues for the cell to evaluate the gradient at the basis functions
         # of the gradient field
         Ferrite.reinit!(cv, cell)
 
         # Now we simply loop over all basis functions of the gradient field and evaluate the gradient
-        for i ∈ 1:num_base_funs
+        for i ∈ 1:Ferrite.getnquadpoints(cv)
             uᵉgradi = Ferrite.function_gradient(cv, i, uᵉ)
             for ds in 1:spatial_dim
                 for df in 1:field_dim

test/runtests.jl

@@ -1,75 +1,179 @@
-using FerriteViz
+using FerriteViz, Ferrite
 using Test
 
-include("../docs/src/ferrite-examples/heat-equation.jl")
+# _test_tolerance(ip::Interpolation{<:Any,<:Any,1}) = 5e-1
+_test_tolerance(ip::Interpolation) = 1e-6 # Float32 computations are involved!
 
-# TODO move this into Ferrite core
+struct MatrixValued <: Ferrite.FieldTrait end
 
-@testset "gradient fields (scalar)" begin
+function Ferrite.function_value(::MatrixValued, fe_v::Ferrite.Values{dim}, q_point::Int, u::AbstractVector) where {dim}
+    n_base_funcs = Ferrite.getn_scalarbasefunctions(fe_v)
+    length(u) == n_base_funcs*dim^2 || Ferrite.throw_incompatible_dof_length(length(u), n_base_funcs)
+    @boundscheck Ferrite.checkquadpoint(fe_v, q_point)
+    val = zero(Tensor{2, dim})
+
+    @inbounds for I ∈ 1:n_base_funcs*dim^2
+         # First flatten to vector
+        i0, c0 = divrem(I - 1, dim^2)
+        i = i0 + 1
+        v = Ferrite.shape_value(fe_v, q_point, i)
+
+        # Then compute matrix index
+        ci0, cj0 = divrem(c0, dim)
+        ci = ci0 + 1
+        cj = cj0 + 1
+
+        val += Ferrite.Tensor{2, dim}((k, l) -> k == ci && l == cj ? v*u[I] : zero(v))
+    end
+
+    return val
+end
+
+@testset "utility operations" begin
     # Check scalar problems
-    for (size, geo, ip) ∈ [(21, Triangle, Lagrange{2,RefTetrahedron,1}()),
-                           (7,Triangle, Lagrange{2,RefTetrahedron,2}()),
-                           (5,Triangle, Lagrange{2,RefTetrahedron,3}()),
-                           #(21,Tetrahedron, Lagrange{3,RefTetrahedron,1}()), #converges rather slowly in the gradient
-                           (7,Tetrahedron, Lagrange{3,RefTetrahedron,2}()),
-                           (21,Quadrilateral, Lagrange{2,RefCube,1}()),
-                           (7,Quadrilateral, Lagrange{2,RefCube,2}()),
-                           #(21,Hexahedron, Lagrange{3,RefCube,1}()), # slows down the pipeline quite a bit, so left out
-                           (7,Hexahedron, Lagrange{3,RefCube,2}())]
-        @testset "($size, $geo, $ip)" begin
-            # Compute solution
-            dh, u = manufactured_heat_problem(geo, ip, size)
-            ip_geo = Ferrite.default_interpolation(typeof(Ferrite.getcells(FerriteViz.getgrid(dh), 1)))
+    for (num_elements_per_dim, geo, ip) ∈ [
+                        #    (4,Triangle, Lagrange{2,RefTetrahedron,1}()),
+                           (2,Triangle, Lagrange{2,RefTetrahedron,2}()),
+                           (2,Triangle, Lagrange{2,RefTetrahedron,3}()),
+                        #    (5,Tetrahedron, Lagrange{3,RefTetrahedron,1}()),
+                           (2,Tetrahedron, Lagrange{3,RefTetrahedron,2}()),
+                        #    (4,Quadrilateral, Lagrange{2,RefCube,1}()),
+                           (2,Quadrilateral, Lagrange{2,RefCube,2}()),
+                        #    (4,Hexahedron, Lagrange{3,RefCube,1}()),
+                           (2,Hexahedron, Lagrange{3,RefCube,2}())
+        ]
+        @testset "scalar($num_elements_per_dim, $geo, $ip)" begin
+            # Get solution
             dim = Ferrite.getdim(ip)
-            qr = QuadratureRule{dim, Ferrite.getrefshape(ip)}(1)
-
-            # Check solution
-            cellvalues = CellScalarValues(qr, Ferrite.getfieldinterpolation(dh, 1), ip_geo);
-            for cell in CellIterator(dh)
-                reinit!(cellvalues, cell)
-                n_basefuncs = getnbasefunctions(cellvalues)
-                coords = getcoordinates(cell)
-                uₑ = u[celldofs(cell)]
-                for q_point in 1:getnquadpoints(cellvalues)
-                    x = spatial_coordinate(cellvalues, q_point, coords)
-                    for i in 1:n_basefuncs
-                        uₐₙₐ    = prod(cos, x*π/2)
-                        uₐₚₚᵣₒₓ = function_value(cellvalues, q_point, uₑ)
-                        @test isapprox(uₐₙₐ, uₐₚₚᵣₒₓ; atol=1e-2)
+            grid = generate_grid(geo, ntuple(x->num_elements_per_dim, dim));
+
+            dh = DofHandler(grid)
+            add!(dh, :u, ip)
+            close!(dh);
+
+            u = Vector{Float64}(undef, ndofs(dh))
+            f_ana(x::Union{Vec{2},FerriteViz.GeometryBasics.Point{2}}) = 0.5x[1]^2 - 2x[2]^2 + x[1]*x[2]
+            f_ana(x::Union{Vec{3},FerriteViz.GeometryBasics.Point{3}}) = -x[1]^2 + 0.3*x[2]^2 + 2*x[3]^2 + 5x[1]*x[2] -   2x[1]*x[3] + 0.1x[3]*x[2]
+            Ferrite.apply_analytical!(u, dh, :u, f_ana)
+
+            @testset "solution fields" begin
+                plotter = FerriteViz.MakiePlotter(dh,u)
+                data = FerriteViz.transfer_solution(plotter,u,process=x->x)[:,1]
+                visible_nodes = .!isnan.(data)# TODO add API
+                @test all(isapprox.(data[visible_nodes], f_ana.(plotter.physical_coords[visible_nodes]); atol=_test_tolerance(ip)))
+            end
+
+            # Compute gradient/flux field
+            @testset "gradient fields" begin
+                (dh_grad, u_grad) = FerriteViz.interpolate_gradient_field(dh, u, :u)
+
+                # Check gradient of solution
+                @testset "interpolate_gradient_field" begin
+                    qr = QuadratureRule{dim,Ferrite.getrefshape(ip)}(2) # TODO sample random point
+                    ip_geo = Ferrite.default_interpolation(geo)
+                    ip_grad = Ferrite.getfieldinterpolation(dh_grad, Ferrite.find_field(dh_grad, :gradient))
+                    cellvalues_grad = Ferrite.CellVectorValues(qr, ip_grad, ip_geo)
+                    for cell in CellIterator(dh_grad)
+                        reinit!(cellvalues_grad, cell)
+                        coords = getcoordinates(cell)
+                        uₑ = @views u_grad[celldofs(cell)]
+                        for q_point in 1:getnquadpoints(cellvalues_grad)
+                            x = spatial_coordinate(cellvalues_grad, q_point, coords)
+                            uₐₚₚᵣₒₓ = function_value(cellvalues_grad, q_point, uₑ)
+                            uₐₙₐ = Tensors.gradient(f_ana, x)
+                            @test all(isapprox.(uₐₙₐ, uₐₚₚᵣₒₓ;atol=_test_tolerance(ip)))
+                        end
+                    end
+                end
+
+                # Check for correct transfer
+                @testset "transfer_solution" begin
+                    plotter_grad = FerriteViz.MakiePlotter(dh_grad,u_grad)
+                    data_grad = FerriteViz.transfer_solution(plotter_grad,u_grad; field_name=:gradient, process=x->x)
+                    visible_nodes_grad = .!isnan.(data_grad)
+                    for i ∈ 1:size(data_grad, 1)
+                        !visible_nodes_grad[i] && continue
+                        x = Vec{dim,Float64}(plotter_grad.physical_coords[i])
+                        ∇uₐₚₚᵣₒₓ = Vec{dim,Float64}(data_grad[i,:])
+                        ∇uₐₙₐ = Tensors.gradient(f_ana, x)
+                        @test all(isapprox.(∇uₐₚₚᵣₒₓ, ∇uₐₙₐ; atol=_test_tolerance(ip)))
                     end
                 end
             end
+        end
+
+        @testset "vector($num_elements_per_dim, $geo, $ip)" begin
+            # Get solution
+            dim = Ferrite.getdim(ip)
+            grid = generate_grid(geo, ntuple(x->num_elements_per_dim, dim));
+
+            dh = DofHandler(grid)
+            add!(dh, :u, dim, ip)
+            close!(dh);
+
+            # Some test functions with rather complicated gradients
+            f_ana(x::Union{Vec{3},FerriteViz.GeometryBasics.Point{3}}) = Vec{3}((
+                -x[1]^2 + 0.3*x[2]^2 + 2*x[3]^2 + 5x[1]*x[2] -   2x[1]*x[3] + 0.1x[3]*x[2],
+                 x[1]^2 - 0.3*x[2]^2 + 1*x[3]^2 - 5x[1]*x[2] +   2x[1]*x[3]               ,
+                          1.3*x[2]^2 - 2*x[3]^2 + 5x[1]*x[2] - 0.7x[1]*x[3] - 0.1x[3]*x[2],
+            ))
+            f_ana(x::Union{Vec{2},FerriteViz.GeometryBasics.Point{2}}) = Vec{2}((
+                -x[1]^2 + 0.3*x[2]^2 +   5x[1]*x[2],
+                 x[1]^2 + 2.3*x[2]^2 - 0.1x[1]*x[2],
+            ))
+            u = Vector{Float64}(undef, ndofs(dh))
+            Ferrite.apply_analytical!(u, dh, :u, f_ana)
+
+            @testset "solution fields" begin
+                plotter = FerriteViz.MakiePlotter(dh,u)
+                data = FerriteViz.transfer_solution(plotter,u,process=x->x)
+                visible_nodes = .!isnan.(data)# TODO add API
+                for i ∈ 1:size(data, 1)
+                    !visible_nodes[i] && continue
+                    uₐₚₚᵣₒₓ = Vec{dim}(data[i,:])
+                    uₐₙₐ = f_ana(Vec{dim}(plotter.physical_coords[i]))
+                    @test all(isapprox.(uₐₚₚᵣₒₓ, uₐₙₐ; atol=_test_tolerance(ip)))
+                end
+            end
 
             # Compute gradient/flux field
-            (dh_grad, u_grad) = FerriteViz.interpolate_gradient_field(dh, u, :u)
-
-            # Check gradient of solution
-            cellvalues_grad = CellVectorValues(qr, Ferrite.getfieldinterpolation(dh_grad, 1), ip_geo);
-            for cell in CellIterator(dh_grad)
-                reinit!(cellvalues_grad, cell)
-                n_basefuncs = getnbasefunctions(cellvalues_grad)
-                coords = getcoordinates(cell)
-                uₑ = u_grad[celldofs(cell)]
-                for q_point in 1:getnquadpoints(cellvalues_grad)
-                    x = spatial_coordinate(cellvalues_grad, q_point, coords)
-                    for i ∈ 1:n_basefuncs
-                        uₐₚₚᵣₒₓ = function_value(cellvalues_grad, q_point, uₑ)
-                        for d ∈ 1:dim
-                            uₐₙₐ = π/2
-                            for j ∈ 1:(d-1)
-                                uₐₙₐ *= cos(x[j]*π/2)
-                            end
-                            uₐₙₐ *= -sin(x[d]*π/2)
-                            for j ∈ (d+1):dim
-                                uₐₙₐ *= cos(x[j]*π/2)
-                            end
-                            @test isapprox(uₐₙₐ, uₐₚₚᵣₒₓ[d]; atol=1e-1)
+            @testset "gradient fields" begin
+                (dh_grad, u_grad) = FerriteViz.interpolate_gradient_field(dh, u, :u)
+
+                # Check gradient of solution
+                @testset "interpolate_gradient_field" begin
+                    qr = QuadratureRule{dim,Ferrite.getrefshape(ip)}(2) # TODO sample random point
+                    ip_geo = Ferrite.default_interpolation(geo)
+                    ip_grad = Ferrite.getfieldinterpolation(dh_grad, Ferrite.find_field(dh_grad, :gradient))
+                    cellvalues_grad = Ferrite.CellScalarValues(qr, ip_grad, ip_geo)
+                    for cell in CellIterator(dh_grad)
+                        reinit!(cellvalues_grad, cell)
+                        coords = getcoordinates(cell)
+                        uₑ = @views u_grad[celldofs(cell)]
+                        for q_point in 1:getnquadpoints(cellvalues_grad)
+                            x = spatial_coordinate(cellvalues_grad, q_point, coords)
+                            ∇uₐₚₚᵣₒₓ = function_value(MatrixValued(), cellvalues_grad, q_point, uₑ)
+                            ∇uₐₙₐ = Tensors.gradient(f_ana, x)
+                            @test all(isapprox.(∇uₐₙₐ, ∇uₐₚₚᵣₒₓ;atol=_test_tolerance(ip)))
                         end
                     end
                 end
+
+                # Check for correct transfer
+                @testset "transfer_solution" begin
+                    plotter_grad = FerriteViz.MakiePlotter(dh_grad,u_grad)
+                    data_grad = FerriteViz.transfer_solution(plotter_grad,u_grad; field_name=:gradient, process=x->x)
+                    visible_nodes_grad = .!isnan.(data_grad)
+                    for i ∈ 1:size(data_grad, 1)
+                        !visible_nodes_grad[i] && continue
+                        x = Vec{dim,Float64}(plotter_grad.physical_coords[i])
+                        # Transpose because constructed from Vector and not from Tuple :)
+                        ∇uₐₚₚᵣₒₓ = transpose(Tensor{2,dim,Float64,2*dim}(data_grad[i,:]))
+                        ∇uₐₙₐ = Tensors.gradient(f_ana, x)
+                        @test all(isapprox.(∇uₐₙₐ, ∇uₐₚₚᵣₒₓ; atol=_test_tolerance(ip)))
+                    end
+                end
             end
         end
     end
 end
-
-# TODO add gradient field test for vector valued problems via manufactured solution