Postprocessing

L2Projector(grid::AbstractGrid)

Initiate an L2Projector for projecting quadrature data onto a function space. To define the function space, add interpolations for differents cell sets with add! before close!ing the projector, see the example below.

The L2Projector acts as the integrated left hand side of the projection equation: Find projection $u \in U_h(\Omega) \subset L_2(\Omega)$ such that

\[\int v u \ \mathrm{d}\Omega = \int v f \ \mathrm{d}\Omega \quad \forall v \in U_h(\Omega),\]

where $f \in L_2(\Omega)$ is the data to project. The function space $U_h(\Omega)$ is the finite element approximation given by the interpolations add!ed to the L2Projector.

Example

proj = L2Projector(grid)
qr_quad = QuadratureRule{RefQuadrilateral}(2)
add!(proj, quad_set, Lagrange{RefQuadrilateral, 1}(); qr_rhs = qr_quad)
qr_tria = QuadratureRule{RefTriangle}(1)
add!(proj, tria_set, Lagrange{RefTriangle, 1}(); qr_rhs = qr_tria)
close!(proj)

vals = Dict{Int, Vector{Float64}}() # Can also be Vector{Vector},
                                    # indexed with cellnr
for (set, qr) in ((quad_set, qr_quad), (tria_set, qr_tria))
    nqp = getnquadpoints(qr)
    for cellnr in set
        vals[cellnr] = rand(nqp)
    end
end

projected = project(proj, vals)

where projected can be used in e.g. evaluate_at_points with the PointEvalHandler, or with evaluate_at_grid_nodes.

add!(proj::L2Projector, set::AbstractVecOrSet{Int}, ip::Interpolation;
    qr_rhs, [qr_lhs])

Add an interpolation ip on the cells in set to the L2Projector proj.

• qr_rhs sets the quadrature rule used to later integrate the right-hand-side of the projection equation, when calling project. It should match the quadrature points used when creating the quadrature-point variables to project.
• The optional qr_lhs sets the quadrature rule used to integrate the left-hand-side of the projection equation, and defaults to a quadrature rule that integrates the mass-matrix exactly for the given interpolation ip.

close!(proj::L2Projector)

Close proj which assembles and calculates the left-hand-side of the projection equation, before doing a Cholesky factorization of the mass-matrix.

L2Projector(ip::Interpolation, grid::AbstractGrid; [qr_lhs], [set])

A quick way to initiate an L2Projector, add an interpolation ip on the set to it, and then close! it so that it can be used to project. The optional keyword argument set defaults to all cells in the grid, while qr_lhs defaults to a quadrature rule that integrates the mass matrix exactly for the interpolation ip.

project(proj::L2Projector, vals, [qr_rhs::QuadratureRule])

Makes a L2 projection of data vals to the nodes of the grid using the projector proj (see L2Projector).

project integrates the right hand side, and solves the projection $u$ from the following projection equation: Find projection $u \in U_h(\Omega) \subset L_2(\Omega)$ such that

\[\int v u \ \mathrm{d}\Omega = \int v f \ \mathrm{d}\Omega \quad \forall v \in U_h(\Omega),\]

where $f \in L_2(\Omega)$ is the data to project. The function space $U_h(\Omega)$ is the finite element approximation given by the interpolations in proj.

The data vals should be an AbstractVector or AbstractDict that is indexed by the cell number. Each index in vals should give an AbstractVector with one element for each cell quadrature point.

If proj was created by calling L2Projector(ip, grid, set), qr_rhs must be given. Otherwise, this is added for each domain when calling add!(proj, args...).

Alternatively, vals can be a matrix, with the column index referring the cell number, and the row index corresponding to quadrature point number. Example (scalar) input data:

vals = [
    [0.44, 0.98, 0.32], # data for quadrature point 1, 2, 3 of element 1
    [0.29, 0.48, 0.55], # data for quadrature point 1, 2, 3 of element 2
    # ...
]

or equivalent in matrix form:

vals = [
    0.44 0.29 # ...
    0.98 0.48 # ...
    0.32 0.55 # ...
]

Supported data types to project are Numbers and AbstractTensors.

evaluate_at_grid_nodes(dh::AbstractDofHandler, u::AbstractVector{T}, fieldname::Symbol) where T

Evaluate the approximated solution for field fieldname at the node coordinates of the grid given the Dof handler dh and the solution vector u.

Return a vector of length getnnodes(grid) where entry i contains the evaluation 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 NaNs.

PointEvalHandler(grid::Grid, points::AbstractVector{Vec{dim,T}}; kwargs...) where {dim, T}

The PointEvalHandler can be used for function evaluation in arbitrary points in the domain – not just in quadrature points or nodes.

The constructor takes a grid and a vector of coordinates for the points. The PointEvalHandler computes i) the corresponding cell, and ii) the (local) coordinate within the cell, for each point. The fields of the PointEvalHandler are:

• cells::Vector{Union{Int,Nothing}}: vector with cell IDs for the points, with nothing for points that could not be found.
• local_coords::Vector{Union{Vec,Nothing}}: vector with the local coordinates (i.e. coordinates in the reference configuration) for the points, with nothing for points that could not be found.

There are two ways to use the PointEvalHandler to evaluate functions:

• evaluate_at_points: can be used when the function is described by i) a dh::DofHandler + uh::Vector (for example the FE-solution), or ii) a p::L2Projector + ph::Vector (for projected data).
• Iteration with PointIterator + PointValues: can be used for more flexible evaluation in the points, for example to compute gradients.

evaluate_at_points(ph::PointEvalHandler, dh::AbstractDofHandler, dof_values::Vector{T}, [fieldname::Symbol]) where T
evaluate_at_points(ph::PointEvalHandler, proj::L2Projector, dof_values::Vector{T}) where T

Return a Vector{T} (for a 1-dimensional field) or a Vector{Vec{fielddim, T}} (for a vector field) with the field values of field fieldname in the points of the PointEvalHandler. The fieldname can be omitted if only one field is stored in dh. The field values are computed based on the dof_values and interpolated to the local coordinates by the function interpolation of the corresponding field stored in the AbstractDofHandler or the L2Projector.

Points that could not be found in the domain when constructing the PointEvalHandler will have NaNs for the corresponding entries in the output vector.

PointValues(cv::CellValues)
PointValues([::Type{T}], func_interpol::Interpolation, [geom_interpol::Interpolation])

Similar to CellValues but with a single updateable "quadrature point". PointValues are used for evaluation of functions/gradients in arbitrary points of the domain together with a PointEvalHandler.

PointValues can be created from CellValues, or from the interpolations directly.

PointValues are reinitialized like other CellValues, but since the local reference coordinate of the "quadrature point" changes this needs to be passed to reinit!, in addition to the element coordinates: reinit!(pv, coords, local_coord). Alternatively, it can be reinitialized with a PointLocation when iterating a PointEvalHandler with a PointIterator.

For function/gradient evaluation, PointValues are used in the same way as CellValues, i.e. by using function_value, function_gradient, etc, with the exception that there is no need to specify the quadrature point index (since PointValues only have 1, this is the default).

PointIterator(ph::PointEvalHandler)

Create an iterator over the points in the PointEvalHandler. The elements of the iterator are either a PointLocation, if the corresponding point could be found in the grid, or nothing, if the point was not found.

A PointLocation can be used to query the cell ID with the cellid function, and can be used to reinitialize PointValues with reinit!.

Examples

ph = PointEvalHandler(grid, points)

for point in PointIterator(ph)
    point === nothing && continue # Skip any points that weren't found
    reinit!(pointvalues, point)   # Update pointvalues
    # ...
end

PointLocation

Element of a PointIterator, typically used to reinitialize PointValues. Fields:

• cid::Int: ID of the cell containing the point
• local_coord::Vec: the local (reference) coordinate of the point
• coords::Vector{Vec}: the coordinates of the cell

VTKGridFile(filename::AbstractString, grid::AbstractGrid; kwargs...)
VTKGridFile(filename::AbstractString, dh::DofHandler; kwargs...)

Create a VTKGridFile that contains an unstructured VTK grid. The keyword arguments are forwarded to WriteVTK.vtk_grid, see Data Formatting Options

This file handler can be used to to write data with

• write_solution
• write_cell_data
• write_projection
• write_node_data.
• Ferrite.write_cellset
• Ferrite.write_nodeset
• Ferrite.write_constraints

It is necessary to call close(::VTKGridFile) to save the data after writing to the file handler. Using the supported do-block does this automatically:

VTKGridFile(filename, grid) do vtk
    write_solution(vtk, dh, u)
    write_cell_data(vtk, celldata)
end

write_solution(vtk::VTKGridFile, dh::AbstractDofHandler, u::Vector, suffix="")

Save the values at the nodes in the degree of freedom vector u to vtk. Each field in dh will be saved separately, and suffix can be used to append to the fieldname.

u can also contain tensorial values, but each entry in u must correspond to a degree of freedom in dh, see write_node_data for details. Use write_node_data directly when exporting values that are already sorted by the nodes in the grid.

write_projection(vtk::VTKGridFile, proj::L2Projector, vals::Vector, name::AbstractString)

Project vals to the grid nodes with proj and save to vtk.

write_cell_data(vtk::VTKGridFile, celldata::AbstractVector, name::String)

Write the celldata that is ordered by the cells in the grid to the vtk file.

write_node_data(vtk::VTKGridFile, nodedata::Vector{Real}, name)
write_node_data(vtk::VTKGridFile, nodedata::Vector{<:AbstractTensor}, name)

Write the nodedata that is ordered by the nodes in the grid to vtk.

When nodedata contains Tensors.Vecs, each component is exported. Two-dimensional vectors are padded with zeros.

When nodedata contains second order tensors, the index order, [11, 22, 33, 23, 13, 12, 32, 31, 21], follows the default Voigt order in Tensors.jl.

write_cellset(vtk, grid::AbstractGrid)
write_cellset(vtk, grid::AbstractGrid, cellset::String)
write_cellset(vtk, grid::AbstractGrid, cellsets::Union{AbstractVector{String},AbstractSet{String})

Write all cell sets in the grid with name according to their keys and celldata 1 if the cell is in the set, and 0 otherwise. It is also possible to only export a single cellset, or multiple cellsets.

write_nodeset(vtk::VTKGridFile, grid::AbstractGrid, nodeset::String)

Write nodal values of 1 for nodes in nodeset, and 0 otherwise

write_constraints(vtk::VTKGridFile, ch::ConstraintHandler)

Saves the dirichlet boundary conditions to a vtkfile. Values will have a 1 where bcs are active and 0 otherwise

write_cell_colors(vtk::VTKGridFile, grid::AbstractGrid, cell_colors, name="coloring")

Write cell colors (see create_coloring) to a VTK file for visualization.

In case of coloring a subset, the cells which are not part of the subset are represented as color 0.