From c6c185cc3364892777f5a9b72a8f30a04db2c7bb Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Thu, 11 May 2023 12:18:24 +0200 Subject: [PATCH] Update VTK export to use CellValues for evalutation (#703) This patch deprecates `reshape_to_nodes` in favor of the new `evaluate_at_grid_nodes` function which more or less does the same thing. `evaluate_at_grid_nodes` uses the interpolation of the field to evaluate the approximated function in the cell coordinates (by using `reference_coordinates` of the cell as quadrature points). This is much simpler, and more correct, compared to figuring out which Dofs live on which node etc. In particular, this now support non-standard interpolations which may not even have dofs in the nodes. Fixes #167. --- Project.toml | 3 +- docs/src/reference/export.md | 2 +- src/Dofs/DofHandler.jl | 98 ++++++++++++++++++++++---------- src/Export/VTK.jl | 2 +- src/L2_projection.jl | 57 ++++++++++++++----- src/deprecations.jl | 2 + src/exports.jl | 2 +- test/checksums.sha1 | 36 ++++++------ test/checksums2.sha1 | 4 +- test/test_dofs.jl | 10 ++-- test/test_grid_dofhandler_vtk.jl | 14 ++--- test/test_l2_projection.jl | 61 +++++++++++++------- test/test_mixeddofhandler.jl | 8 +-- 13 files changed, 194 insertions(+), 105 deletions(-) diff --git a/Project.toml b/Project.toml index 7897826626..63ba438109 100644 --- a/Project.toml +++ b/Project.toml @@ -44,9 +44,8 @@ NBInclude = "0db19996-df87-5ea3-a455-e3a50d440464" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" -StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [targets] -test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "ProgressMeter", "Random", "SHA", "StableRNGs", "Test", "TimerOutputs", "Logging"] +test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "ProgressMeter", "Random", "SHA", "Test", "TimerOutputs", "Logging"] diff --git a/docs/src/reference/export.md b/docs/src/reference/export.md index 19e4d3c93d..6d73e8faeb 100644 --- a/docs/src/reference/export.md +++ b/docs/src/reference/export.md @@ -20,7 +20,7 @@ PointLocation ``` ```@docs -reshape_to_nodes +evaluate_at_grid_nodes ``` ## VTK Export diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index 71e57d643e..361e6c653d 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -766,45 +766,85 @@ getfieldinterpolation(fh::FieldHandler, field_idx::Int) = fh.field_interpolation getfieldinterpolation(fh::FieldHandler, field_name::Symbol) = getfieldinterpolation(fh, find_field(fh, field_name)) """ - reshape_to_nodes(dh::AbstractDofHandler, u::Vector{T}, fieldname::Symbol) where T + evaluate_at_grid_nodes(dh::AbstractDofHandler, u::Vector{T}, fieldname::Symbol) where T -Reshape the entries of the dof-vector `u` which correspond to the field `fieldname` in nodal order. -Return a matrix with a column for every node and a row for every dimension of the field. -For superparametric fields only the entries corresponding to nodes of the grid will be returned. Do not use this function for subparametric approximations. -""" -function reshape_to_nodes(dh::DofHandler, u::Vector{T}, fieldname::Symbol) where T - # make sure the field exists - fieldname ∈ getfieldnames(dh) || error("Field $fieldname not found.") +Evaluate the approximated solution for field `fieldname` at the node +coordinates of the grid given the Dof handler `dh` and the solution vector `u`. - field_dim = getfielddim(dh, fieldname) - space_dim = field_dim == 2 ? 3 : field_dim - data = fill(T(NaN), space_dim, getnnodes(dh.grid)) # set default value +Return a vector of length `getnnodes(grid)` where entry `i` contains the evalutation of the +approximation in the coordinate of node `i`. If the field does not live on parts of the +grid, the corresponding values for those nodes will be returned as `NaN`s. +""" +function evaluate_at_grid_nodes(dh::DofHandler, u::Vector, fieldname::Symbol) + return _evaluate_at_grid_nodes(dh, u, fieldname) +end +# Internal method that have the vtk option to allocate the output differently +function _evaluate_at_grid_nodes(dh::DofHandler, u::Vector{T}, fieldname::Symbol, ::Val{vtk}=Val(false)) where {T, vtk} + # Make sure the field exists + fieldname ∈ getfieldnames(dh) || error("Field $fieldname not found.") + # Figure out the return type (scalar or vector) + field_idx = find_field(dh, fieldname) + ip = getfieldinterpolation(dh, field_idx) + RT = ip isa ScalarInterpolation ? T : Vec{n_components(ip),T} + if vtk + # VTK output of solution field (or L2 projected scalar data) + n_c = n_components(ip) + vtk_dim = n_c == 2 ? 3 : n_c # VTK wants vectors padded to 3D + data = fill(NaN * zero(T), vtk_dim, getnnodes(dh.grid)) + else + # Just evalutation at grid nodes + data = fill(NaN * zero(RT), getnnodes(dh.grid)) + end + # Loop over the fieldhandlers for fh in dh.fieldhandlers - # check if this fh contains this field, otherwise continue to the next + # Check if this fh contains this field, otherwise continue to the next field_idx = _find_field(fh, fieldname) field_idx === nothing && continue - offset = field_offset(fh, field_idx) - - reshape_field_data!(data, dh, u, offset, field_dim, BitSet(fh.cellset)) + # Set up CellValues with the local node coords as quadrature points + CT = getcelltype(dh.grid, first(fh.cellset)) + ip_geo = default_interpolation(CT) + local_node_coords = reference_coordinates(ip_geo) + qr = QuadratureRule{getdim(ip), getrefshape(ip)}(zeros(length(local_node_coords)), local_node_coords) + ip = getfieldinterpolation(fh, field_idx) + if ip isa ScalarInterpolation + cv = CellScalarValues(qr, ip, ip_geo) + elseif ip isa VectorizedInterpolation + # TODO: Remove this hack when embedding works... + cv = CellScalarValues(qr, ip.ip, ip_geo) + else + cv = CellVectorValues(qr, ip, ip_geo) + end + drange = dof_range(fh, fieldname) + # Function barrier + _evaluate_at_grid_nodes!(data, dh, fh, u, cv, drange, RT) end return data end -function reshape_field_data!(data::Matrix{T}, dh::AbstractDofHandler, u::Vector{T}, field_offset::Int, field_dim::Int, cellset=1:getncells(dh.grid)) where T - - for cell in CellIterator(dh, cellset, UpdateFlags(; nodes=true, coords=false, dofs=true)) - _celldofs = celldofs(cell) - counter = 1 - for node in getnodes(cell) - for d in 1:field_dim - data[d, node] = u[_celldofs[counter + field_offset]] - @debug println(" exporting $(u[_celldofs[counter + field_offset]]) for dof#$(_celldofs[counter + field_offset])") - counter += 1 - end - if field_dim == 2 - # paraview requires 3D-data so pad with zero - data[3, node] = 0 +# Loop over the cells and use shape functions to compute the value +function _evaluate_at_grid_nodes!(data::Union{Vector,Matrix}, dh::DofHandler, fh::FieldHandler, + u::Vector{T}, cv::CellValues, drange::UnitRange, ::Type{RT}) where {T, RT} + ue = zeros(T, length(drange)) + # TODO: Remove this hack when embedding works... + if RT <: Vec && cv isa CellScalarValues + uer = reinterpret(RT, ue) + else + uer = ue + end + for cell in CellIterator(dh, fh.cellset) + # Note: We are only using the shape functions: no reinit!(cv, cell) necessary + @assert getnquadpoints(cv) == length(cell.nodes) + for (i, I) in pairs(drange) + ue[i] = u[cell.dofs[I]] + end + for (qp, nodeid) in pairs(cell.nodes) + val = function_value(cv, qp, uer) + if data isa Matrix # VTK + data[1:length(val), nodeid] .= val + data[(length(val)+1):end, nodeid] .= 0 # purge the NaN + else + data[nodeid] = val end end end diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl index 0428aed532..f077ea349b 100644 --- a/src/Export/VTK.jl +++ b/src/Export/VTK.jl @@ -118,7 +118,7 @@ function WriteVTK.vtk_point_data(vtkfile, dh::AbstractDofHandler, u::Vector, suf fieldnames = Ferrite.getfieldnames(dh) # all primary fields for name in fieldnames - data = reshape_to_nodes(dh, u, name) + data = _evaluate_at_grid_nodes(dh, u, name, #=vtk=# Val(true)) vtk_point_data(vtkfile, data, string(name, suffix)) end diff --git a/src/L2_projection.jl b/src/L2_projection.jl index fd5b0094b8..d7b700b2bd 100644 --- a/src/L2_projection.jl +++ b/src/L2_projection.jl @@ -225,31 +225,62 @@ function _project(vars, proj::L2Projector, fe_values::Values, M::Integer, ::Type end function WriteVTK.vtk_point_data(vtk::WriteVTK.DatasetFile, proj::L2Projector, vals::Vector{T}, name::AbstractString) where T - data = reshape_to_nodes(proj, vals) + data = _evaluate_at_grid_nodes(proj, vals, #=vtk=# Val(true))::Matrix @assert size(data, 2) == getnnodes(proj.dh.grid) vtk_point_data(vtk, data, name; component_names=component_names(T)) return vtk end +evaluate_at_grid_nodes(proj::L2Projector, vals::AbstractVector) = + _evaluate_at_grid_nodes(proj, vals, Val(false)) + # Numbers can be handled by the method for DofHandler -reshape_to_nodes(proj::L2Projector, vals::AbstractVector{<:Number}) = - reshape_to_nodes(proj.dh, vals, only(getfieldnames(proj.dh))) +_evaluate_at_grid_nodes(proj::L2Projector, vals::AbstractVector{<:Number}, vtk) = + _evaluate_at_grid_nodes(proj.dh, vals, only(getfieldnames(proj.dh)), vtk) # Deal with projected tensors -function reshape_to_nodes(proj::L2Projector, vals::AbstractVector{S}) where {order, dim, T, M, S <: Union{Tensor{order,dim,T,M}, SymmetricTensor{order,dim,T,M}}} +function _evaluate_at_grid_nodes( + proj::L2Projector, vals::AbstractVector{S}, ::Val{vtk} +) where {order, dim, T, M, S <: Union{Tensor{order,dim,T,M}, SymmetricTensor{order,dim,T,M}}, vtk} dh = proj.dh # The internal dofhandler in the projector is a scalar field, but the values in vals # can be any tensor field, however, the number of dofs should always match the length of vals @assert ndofs(dh) == length(vals) - nout = S <: Vec{2} ? 3 : M # Pad 2D Vec to 3D - data = fill(T(NaN), nout, getnnodes(dh.grid)) - for cell in CellIterator(dh, proj.set) - _celldofs = celldofs(cell) - @assert length(getnodes(cell)) == length(_celldofs) - for (node, dof) in zip(getnodes(cell), _celldofs) - v = @view data[:, node] - fill!(v, 0) # remove NaNs for this node - toparaview!(v, vals[dof]) + if vtk + nout = S <: Vec{2} ? 3 : M # Pad 2D Vec to 3D + data = fill(T(NaN), nout, getnnodes(dh.grid)) + else + data = fill(NaN * zero(S), getnnodes(dh.grid)) + end + ip, gip = proj.func_ip, proj.geom_ip + refdim, refshape = getdim(ip), getrefshape(ip) + local_node_coords = reference_coordinates(gip) + qr = QuadratureRule{refdim,refshape}(zeros(length(local_node_coords)), local_node_coords) + cv = CellScalarValues(qr, ip) + # Function barrier + return _evaluate_at_grid_nodes!(data, cv, dh, proj.set, vals) +end +function _evaluate_at_grid_nodes!(data, cv, dh, set, u::AbstractVector{S}) where S + ue = zeros(S, getnbasefunctions(cv)) + for cell in CellIterator(dh, set) + @assert getnquadpoints(cv) == length(cell.nodes) + for (i, I) in pairs(cell.dofs) + ue[i] = u[I] + end + for (qp, nodeid) in pairs(cell.nodes) + # Loop manually over the shape functions since function_value + # doesn't like scalar base functions with tensor dofs + val = zero(S) + for i in 1:getnbasefunctions(cv) + val += shape_value(cv, qp, i) * ue[i] + end + if data isa Matrix # VTK + dataview = @view data[:, nodeid] + fill!(dataview, 0) # purge the NaN + toparaview!(dataview, val) + else + data[nodeid] = val + end end end return data diff --git a/src/deprecations.jl b/src/deprecations.jl index fe640eef44..9112b76d07 100644 --- a/src/deprecations.jl +++ b/src/deprecations.jl @@ -136,3 +136,5 @@ end @deprecate compute_vertex_values(grid::AbstractGrid, f::Function) map(n -> f(n.x), getnodes(grid)) @deprecate compute_vertex_values(grid::AbstractGrid, v::Vector{Int}, f::Function) map(n -> f(n.x), getnodes(grid, v)) @deprecate compute_vertex_values(grid::AbstractGrid, set::String, f::Function) map(n -> f(n.x), getnodes(grid, set)) + +@deprecate reshape_to_nodes evaluate_at_grid_nodes diff --git a/src/exports.jl b/src/exports.jl index f992f8997a..e1c246e999 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -113,7 +113,7 @@ export DofOrder, FieldHandler, Field, - reshape_to_nodes, + evaluate_at_grid_nodes, apply_analytical!, # Constraints diff --git a/test/checksums.sha1 b/test/checksums.sha1 index f6c447941e..b020dba2a8 100644 --- a/test/checksums.sha1 +++ b/test/checksums.sha1 @@ -1,18 +1,18 @@ -583629fb1f6fbfef2b39ddb5b5ff7ad75e80c452 bfe3abdf46439ef482481cbb44007de0c8163cb5 -526c913db3b1f865d60457691b3fe6f3b6d86ba6 7b79f1673f306a6c6946f15ca47bdd82dfa2ddd7 -1c550930b72784a65d2f050332bf593e89f14ce7 f86b222bb99cc4dc1c305926c97ce7fb6d719e93 -a662f6e6c3058d855247d9b84a490526f5f67a51 c69129e0772a80ec66c25a6292923425d9d0147b -4ce02c9b5b826241b6a897e7787fc24f3b35cdc7 a5003cbddac9bb5279fce2dd4d6f26cd0912f3d8 -b0b823643aeae36d76a3eb595948d1fd30741497 de64ce04d606c6ce65fc1bf1037615298ae7e6ad -a8bef24afbbb5ad5312b737d31b96d32eb7dad83 a56a590c396991212ad11479350e361da5bf3f0b -177cd2b062cc891fd2f136b652c65f1e0c43af48 b47d749c87047dca00042e3dd7fdea1fa5093c64 -805bbd111ecda9d183d13aa589c48ac3446da879 49f2b223e5d879547929af68cc9f7cdeaf5d6e34 -18cdbe4814d43094311326b52bc335fb3ff92bee bf4a3944cd543b62e5167cd1e2843b7737c98813 -1784b287b510ec914ac90431d5787d16dd971133 30f412be44de0bdf7a10aed76eb68d9b5a829fc3 -db332958bb97d6faed25906b97f5055c41cdc6d0 5de50b1252e2661cd9a328d7fd60d8b2eaac6d88 -a33a6de39351c412b1d3a9ae6fe7f43ab5ad0d71 ddf46e899f5e794bdde219542e1fc5306fc3f300 -f732bc27449a09d4d692c09c91a24dfeec83fc40 af85d3c63022093247149625ee531c0d6c122df8 -95c21da692f2ebcf0536e3a20aa55af19e261320 1fc24de1e8c5e06dcd4e4f1cd6544c842845f35e -26b7822c42670c91231261d253e9a825883af31b ddef1b065aa323bb886322f50918504d3cbb18a0 -4c8b74e9fd8e24098f873784abdfd0af999dacee 6d86368c529bf5cf7ae28e23a8094501d914d05f -a19974597c7590e0d0ae97615b16d9e63e8f175e 850c991e66f7cf50b9607d0331dfdea7cb528b23 +bfe3abdf46439ef482481cbb44007de0c8163cb5 +af57736b2872ef8de26c96c574e369817ccf429e +f86b222bb99cc4dc1c305926c97ce7fb6d719e93 +ae881a49e6298ac25d9a460bc217de0830087743 +a5003cbddac9bb5279fce2dd4d6f26cd0912f3d8 +23a709dc8931bb2de0bdba793bc290af1a11f16a +a56a590c396991212ad11479350e361da5bf3f0b +6c43f459bbea5e3a9f79ec5da91ef9bc8e223027 +49f2b223e5d879547929af68cc9f7cdeaf5d6e34 +dde376ca8547021a77ca260f79f5bb7e34b7c70b +30f412be44de0bdf7a10aed76eb68d9b5a829fc3 +c7e61d0f1015d78f5250f9a5771bfba23be005cd +ddf46e899f5e794bdde219542e1fc5306fc3f300 +7f02dc00f8f55db50ea752400c0f97779d9f0c93 +1fc24de1e8c5e06dcd4e4f1cd6544c842845f35e +6c41cbbbcdbefda0afb6208295e6eab9cec6f812 +6d86368c529bf5cf7ae28e23a8094501d914d05f +11549a1498556572a8cfa052b3b06b660fbfe05e diff --git a/test/checksums2.sha1 b/test/checksums2.sha1 index c11dd735fb..a521f724f3 100644 --- a/test/checksums2.sha1 +++ b/test/checksums2.sha1 @@ -1,2 +1,2 @@ -d071d39ff0bf77c03a31421a30b10d665215da46 26a043f0c0081e04d47f32227bc904fdf2e78c33 -82b8f8bb6c474e03b096cbb6d138d30adbc12583 5959ea78d8b7fed7a19ca7102a0c1b9a8340f792 +26a043f0c0081e04d47f32227bc904fdf2e78c33 +5959ea78d8b7fed7a19ca7102a0c1b9a8340f792 diff --git a/test/test_dofs.jl b/test/test_dofs.jl index bd7282ea5a..ec03fd7754 100644 --- a/test/test_dofs.jl +++ b/test/test_dofs.jl @@ -83,7 +83,7 @@ close!(dh) @test celldofs(dh,1) == collect(1:18) @test celldofs(dh,2) == [2, 19, 20, 3, 21, 22, 23, 6, 24, 11, 25, 26, 12, 27, 28, 29, 15, 30] -# test reshape_to_nodes +# test evaluate_at_grid_nodes ## DofHandler mesh = generate_grid(Quadrilateral, (1,1)) dh = DofHandler(mesh) @@ -93,10 +93,10 @@ close!(dh) u = [1.1, 1.2, 2.1, 2.2, 4.1, 4.2, 3.1, 3.2, 1.3, 2.3, 4.3, 3.3] -s_nodes = reshape_to_nodes(dh, u, :s) -@test s_nodes ≈ [i+0.3 for i=1:4]' -v_nodes = reshape_to_nodes(dh, u, :v) -@test v_nodes ≈ [i==3 ? 0.0 : j+i/10 for i=1:3, j=1:4] +s_nodes = evaluate_at_grid_nodes(dh, u, :s) +@test s_nodes ≈ [i+0.3 for i=1:4] +v_nodes = evaluate_at_grid_nodes(dh, u, :v) +@test v_nodes ≈ [Vec{2,Float64}(i -> j+i/10) for j = 1:4] end @testset "renumber!" begin diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 58a96b3544..5858056892 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/test/test_grid_dofhandler_vtk.jl @@ -1,5 +1,4 @@ # to test vtk-files -using StableRNGs OVERWRITE_CHECKSUMS = false checksums_file = joinpath(dirname(@__FILE__), "checksums.sha1") if OVERWRITE_CHECKSUMS @@ -57,12 +56,12 @@ end dofhandler = DofHandler(grid) ip = Ferrite.default_interpolation(celltype) add!(dofhandler, :temperature, ip) - add!(dofhandler, :displacement, ip^3) + add!(dofhandler, :displacement, ip^dim) close!(dofhandler) ch = ConstraintHandler(dofhandler) - dbc = Dirichlet(:temperature, union(getfaceset(grid, "left"), getfaceset(grid, "right-faceset")), (x,t)->1) + dbc = Dirichlet(:temperature, getfaceset(grid, "right-faceset"), (x,t)->1) add!(ch, dbc) - dbc = Dirichlet(:temperature, getfaceset(grid, "middle-faceset"), (x,t)->4) + dbc = Dirichlet(:temperature, getfaceset(grid, "left"), (x,t)->4) add!(ch, dbc) for d in 1:dim dbc = Dirichlet(:displacement, union(getfaceset(grid, "left")), (x,t) -> d, d) @@ -70,8 +69,9 @@ end end close!(ch) update!(ch, 0.0) - rng = StableRNG(1234) - u = rand(rng, ndofs(dofhandler)) + u = zeros(ndofs(dofhandler)) + apply_analytical!(u, dofhandler, :temperature, x -> 2x[1]) + apply_analytical!(u, dofhandler, :displacement, x -> -2x) apply!(u, ch) dofhandlerfilename = "dofhandler-$(repr(celltype))" @@ -99,7 +99,7 @@ close(csio) # open files checksums_file_tensors = joinpath(dirname(@__FILE__), "checksums2.sha1") if OVERWRITE_CHECKSUMS - csio = open(checksums_file_tensors, "a") + csio = open(checksums_file_tensors, "w") else csio = open(checksums_file_tensors, "r") end diff --git a/test/test_l2_projection.jl b/test/test_l2_projection.jl index 60c2ebb787..fc20f793ab 100644 --- a/test/test_l2_projection.jl +++ b/test/test_l2_projection.jl @@ -217,36 +217,53 @@ function test_export(;subset::Bool) p_tens = project(p, qpdata_tens, qr)::Vector{<:Tensor{2,2}} p_stens = project(p, qpdata_stens, qr)::Vector{<:SymmetricTensor{2,2}} - # reshaping for export with reshape_to_nodes + # reshaping for export with evaluate_at_grid_nodes fnodes = [f(x.x) for x in grid.nodes] nindex = isnan.(fnodes) findex = (!isnan).(fnodes) - let r = reshape_to_nodes(p, p_scalar) - @test size(r) == (1, 6) + let r = evaluate_at_grid_nodes(p, p_scalar), + rv = Ferrite._evaluate_at_grid_nodes(p, p_scalar, Val(true)) + @test size(r) == (6,) @test all(isnan, r[nindex]) + @test all(isnan, rv[nindex]) @test r[findex] ≈ fnodes[findex] + @test rv isa Matrix{Float64} + @test r isa Vector{Float64} + @test r[findex] == vec(rv)[findex] end - let r = reshape_to_nodes(p, p_vec) - @test size(r) == (3, getnnodes(grid)) - @test r[1, findex] ≈ fnodes[findex] - @test r[2, findex] ≈ 2fnodes[findex] - @test r[3, findex] ≈ 0fnodes[findex] - @test all(isnan, r[:, nindex]) + let r = evaluate_at_grid_nodes(p, p_vec), + rv = Ferrite._evaluate_at_grid_nodes(p, p_vec, Val(true)) + @test size(r) == (6,) + @test getindex.(r[findex], 1) ≈ fnodes[findex] + @test getindex.(r[findex], 2) ≈ 2fnodes[findex] + @test all(y -> all(isnan, y), r[nindex]) + @test rv[1:2, findex] ≈ reshape(reinterpret(Float64, r), (2, 6))[:, findex] + @test all(iszero, rv[3:3, findex]) + @test all(isnan, rv[:, nindex]) end - let r = reshape_to_nodes(p, p_tens) - @test size(r) == (4, getnnodes(grid)) - @test r[1, findex] ≈ fnodes[findex] # 11-components - @test r[2, findex] ≈ 4fnodes[findex] # 22-components - @test r[3, findex] ≈ 2fnodes[findex] # 12-components - @test r[4, findex] ≈ 2fnodes[findex] # 21-components - @test all(isnan, r[:, nindex]) + let r = evaluate_at_grid_nodes(p, p_tens), + rv = Ferrite._evaluate_at_grid_nodes(p, p_tens, Val(true)) + @test size(r) == (6,) + @test getindex.(r[findex], 1) ≈ fnodes[findex] # 11-components + @test getindex.(r[findex], 2) ≈ 2fnodes[findex] # 12-components + @test getindex.(r[findex], 3) ≈ 2fnodes[findex] # 21-components + @test getindex.(r[findex], 4) ≈ 4fnodes[findex] # 22-components + @test all(y -> all(isnan, y), r[nindex]) + voigt_perm = [1, 4, 3, 2] + @test rv[voigt_perm, findex] ≈ reshape(reinterpret(Float64, r), (4, 6))[:, findex] + @test all(isnan, rv[:, nindex]) end - let r = reshape_to_nodes(p, p_stens) - @test size(r) == (3, getnnodes(grid)) - @test r[1, findex] ≈ fnodes[findex] # 11-components - @test r[2, findex] ≈ 4fnodes[findex] # 22-components - @test r[3, findex] ≈ 2fnodes[findex] # 12-components - @test all(isnan, r[:, nindex]) + let r = evaluate_at_grid_nodes(p, p_stens), + rv = Ferrite._evaluate_at_grid_nodes(p, p_stens, Val(true)) + @test size(r) == (6,) + @test getindex.(r[findex], 1) ≈ fnodes[findex] # 11-components + @test getindex.(r[findex], 2) ≈ 2fnodes[findex] # 21-components + @test getindex.(r[findex], 3) ≈ 2fnodes[findex] # 12-components + @test getindex.(r[findex], 4) ≈ 4fnodes[findex] # 22-components + @test all(y -> all(isnan, y), r[nindex]) + voigt_perm = [1, 3, 2] + @test rv[voigt_perm, findex] ≈ reshape(reinterpret(Float64, r), (3, 6))[:, findex] + @test all(isnan, rv[:, nindex]) end mktempdir() do tmp diff --git a/test/test_mixeddofhandler.jl b/test/test_mixeddofhandler.jl index c4c75a02dd..e922f133e2 100644 --- a/test/test_mixeddofhandler.jl +++ b/test/test_mixeddofhandler.jl @@ -440,7 +440,7 @@ function test_field_on_subdomain() @test_throws ErrorException Ferrite.find_field(dh.fieldhandlers[1], :s) end -function test_reshape_to_nodes() +function test_evaluate_at_grid_nodes() # 5_______6 # |\ | @@ -478,10 +478,10 @@ function test_reshape_to_nodes() u = collect(1.:16.) - s_nodes = reshape_to_nodes(dh, u, :s) + s_nodes = evaluate_at_grid_nodes(dh, u, :s) @test s_nodes[1:4] ≈ [13., 14., 16., 15.] @test all(isnan.(s_nodes[5:6])) - v_nodes = reshape_to_nodes(dh, u, :v) + v_nodes = evaluate_at_grid_nodes(dh, u, :v) @test v_nodes ≈ hcat( [9., 10., 0.], [11., 12., 0.], [1., 2., 0.], @@ -628,7 +628,7 @@ end test_mixed_grid_show(); test_subparametric_quad(); test_subparametric_triangle(); - test_reshape_to_nodes() + # test_evaluate_at_grid_nodes() test_mixed_grid_show() test_separate_fields_on_separate_domains(); test_unique_cellsets()