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.
Ferrite-FEM · May 11, 2023 · c6c185c · c6c185c
1 parent ade9292
commit c6c185c
Show file tree

Hide file tree

Showing 13 changed files with 194 additions and 105 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -1,2 +1,2 @@
-d071d39ff0bf77c03a31421a30b10d665215da46 26a043f0c0081e04d47f32227bc904fdf2e78c33
-82b8f8bb6c474e03b096cbb6d138d30adbc12583 5959ea78d8b7fed7a19ca7102a0c1b9a8340f792
+26a043f0c0081e04d47f32227bc904fdf2e78c33
+5959ea78d8b7fed7a19ca7102a0c1b9a8340f792
diff --git 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
@@ -1,5 +1,4 @@
 # to test vtk-files
-using StableRNGs
 OVERWRITE_CHECKSUMS = false
 checksums_file = joinpath(dirname(@__FILE__), "checksums.sha1")
 if OVERWRITE_CHECKSUMS
@@ -57,21 +56,22 @@ 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)
             add!(ch, dbc)
         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