Example integrating with DelaunayTriangulation.jl

[fredrikekre]

It would be nice to have an example where https://github.com/JuliaGeometry/DelaunayTriangulation.jl is used for the grid generation. Perhaps we can spice up one of the existing examples that uses a hypercube grid to something more fun.

The author of DelaunayTriangulation.jl (@DanielVandH) have offered to help out with any questions that may arise.

[DanielVandH]

I could probably have a go at this myself with a bit of help understanding the Grid struct. I may be the wrong person for it since I've never actually used this package, but maybe it's useful. A realistic example could be a curve-bounded domain, e.g. https://juliageometry.github.io/DelaunayTriangulation.jl/stable/tutorials/curve_bounded/#A-Complicated-Multiply-Connected-Disjoint-Domain (not this domain exactly since I don't imagine it's very practical, but just to show what could be done).

For example of what I'm having difficulty understanding, with

tri = triangulate_rectangle(0, 1, 0, 1, 4, 4, single_boundary = false) # [0, 1] × [0, 1], 20 points in each direction
triplot(tri)

each_solid_triangle(tri) |> collect

18-element Vector{Tuple{Int64, Int64, Int64}}:
 (9, 10, 13)
 (1, 2, 5)
 (7, 4, 8)
 (6, 3, 7)
 (7, 8, 11)
 (11, 8, 12)
 (9, 6, 10)
 (2, 3, 6)
 (10, 7, 11)
 (11, 12, 15)
 (3, 4, 7)
 (5, 2, 6)
 (6, 7, 10)
 (13, 10, 14)
 (5, 6, 9)
 (14, 11, 15)
 (10, 11, 14)
 (15, 12, 16)

The cells can be easily obtained by iterating over each_solid_triangle(tri), and the nodes from each_solid_vertex and get_point. (If you are unfamiliar, the corresponding _ghost iterators are used for boundary information, as used in e.g. FiniteVolumeMethod.jl.) I don't really understand what (:cellsets, :nodesets, :facetsets, :vertexsets) are defined as exactly, nor their ordering (which I imagine is important given the use of OrderedSets). For example,

grid = generate_grid(Triangle, (4, 4))

julia> grid.cellsets
Dict{String, OrderedCollections.OrderedSet{Int64}}()

julia> grid.nodesets
Dict{String, OrderedCollections.OrderedSet{Int64}}()

julia> grid.facetsets
Dict{String, OrderedCollections.OrderedSet{FacetIndex}} with 4 entries:
  "left"   => OrderedSet{FacetIndex}([FacetIndex((1, 3)), FacetIndex((9, 3)), FacetIndex((17, 3)), FacetIndex((25, 3))])
  "bottom" => OrderedSet{FacetIndex}([FacetIndex((1, 1)), FacetIndex((3, 1)), FacetIndex((5, 1)), FacetIndex((7, 1))])
  "right"  => OrderedSet{FacetIndex}([FacetIndex((8, 1)), FacetIndex((16, 1)), FacetIndex((24, 1)), FacetIndex((32, 1))])
  "top"    => OrderedSet{FacetIndex}([FacetIndex((26, 2)), FacetIndex((28, 2)), FacetIndex((30, 2)), FacetIndex((32, 2))])

julia> grid.vertexsets
Dict{String, OrderedCollections.OrderedSet{VertexIndex}}()

I'm not sure what I'm actually looking at here. I did see https://ferrite-fem.github.io/Ferrite.jl/dev/topics/grid/, but I don't really follow it unfortunately.

I'm also assuming that the aim would be to put all the information into Grid at least, which is probably simplest just for a one-off example. Another way would be to define a wrapper

struct DelaunayGrid{T<:Triangulation} <: AbstractGrid
    triangulation::T
end

and use that as the interface. Could be for a useful package extension in the future maybe.

[KnutAM]

Cool!

I don't really understand what (:cellsets, :nodesets, :facetsets, :vertexsets) are defined as exactly

These sets define subsets of the grid, to allow performing operations on only the subset.
The two most common to be filled are cellsets and facetsets. cellsets simply gives a list of cell numbers that can be used to perform specific operations on part of the domain (for example having different physical properties). For making this example, I don't think that will be necessary.

The facetsets are typically used to define boundary conditions for the simulation. The FacetIndex is basically a tuple of (cell_number, local_facet_number), where in 2d the local_facet_number implies which of the edges of the cell the FacetIndex refers to.

However, you do not have to populate this either, since that can be done by calling add_facetset! or add_boundary_facetset! after having created the grid. With your boundary iterator it might be possible to make more efficient, but that could be added later!

The ordering of any of these should have no effect on the simulation results, but may affect performance during assembly and the resulting finite element matrix structure (for cellset).

[DanielVandH]

The FacetIndex is basically a tuple of (cell_number, local_facet_number), where in 2d the local_facet_number implies which of the edges of the cell the FacetIndex refers to.

I see. In grid.facetsets in my example, I see that the left boundary has a FacetIndex(1, 3), so

julia> cn, lcn = FacetIndex(1, 3)
FacetIndex((1, 3))

julia> tri = grid.cells[cn]
Triangle((1, 2, 6))

julia> i, j, k = tri.nodes;

julia> ((i, j), (j, k), (k, i))[lcn]
(6, 1)

julia> grid.nodes[[6, 1]]
2-element Vector{Node{2, Float64}}:
 Node{2, Float64}([-1.0, -0.5])
 Node{2, Float64}([-1.0, -1.0])

shows that the associated triangle is (1, 2, 6) and the edge (6, 1) belongs to the left boundary. I'll try and write a simple example with this understanding.

[termi-official]

Thanks for the input @DanielVandH ! Let me comment on this, trying to help you understand decisions made on our side.

Currently we essentially assume that we can query the following information from the grid efficiently:

1. What is element i and which node indices correspond to the element, which is stored in grid.cells
2. Which coordinates is associated with the node grid.nodes
3. Which subdomains does the grid have (think e.g. about where is the fluid, where is the solid and where is the inflow) which is stored in grid.cellsets and so on

For subdomain information we use ordered grids to utilize caches as efficient as possible. If you use a Set of indices, then the iteration order does not follow the natural iteration order of your cells. In benchmarks we have shown that the iteration over the OrderedSet gives you quite a bit of speed up (4x and more, depending on how hard you try to optimize cache utilization).

For the indices itself we have the cells and nodes, which are just integers (as they index into grid.cells and grid.nodes directly) and we have furthermore indices for subsets of cells, e.g. facets. Since our grid does not materialize facets and ridges explicitly, we have for now settled on this trick to still be able to index into these entities.

If there are further questions do not hesitate to ask.

[KnutAM]

i, j, k = tri.nodes
((i, j), (j, k), (k, i))[lcn]

We have the function

Ferrite.facets(tri)[lcn]

which does this operation btw.
As long as you only have linear triangles, this will return all nodes on the facet, as this corresponds to the vertices of the triangle. For 2nd order triangles, there is a difference between nodes (all points defining the triangle) and the vertices (only the corners of the triangles).

[DanielVandH]

Thanks for your comments. Here are some functions that are useful for this. I might be using some internals of DelaunayTriangulation here, I'd have to double check. If I am, I'll make them public.

First, getting the cells and nodes is very easy. For get_cells, I also a return a mapping that takes a triangle to its position in cells, since my methods for getting the facetsets would be pretty slow otherwise.

const DT = DelaunayTriangulation
function get_cells(tri::Triangulation)
    cells = Vector{Triangle}(undef, num_solid_triangles(tri))
    cell_map = Dict{NTuple{3,Int},Int}()
    for (cell_number, triangle) in enumerate(each_solid_triangle(tri))
        ijk = triangle_vertices(triangle)
        uvw = DT.sort_triangle(ijk) # It's useful to know the order of the vertices
        cells[cell_number] = Triangle(uvw)
        cell_map[uvw] = cell_number
    end
    return cells, cell_map
end
function get_nodes(tri::Triangulation)
    nodes = Vector{Node{2,Float64}}(undef, num_solid_vertices(tri))
    for (node_number, vertex) in enumerate(each_point_index(tri)) # was originally each_solid_vertex, but that breaks if there are any points not in the triangulation.
        p = get_point(tri, vertex)
        x, y = getxy(p)
        nodes[node_number] = Node((x, y))
    end
    return nodes
end

For getting the boundary information, there are actually several ways to do this inside DelaunayTriangulation.jl.

• Iterate over each_boundary_edge. This is unordered.
• Use each_ghost_vertex together with get_adjacent2vertex. This is also unordered, but each of the facet sets can be built one at a time.
• If we wanted to iterate over the edges so that order is respected, we can use get_boundary_nodes. This requires a bit of care to get correct.

I think the third is probably fastest to work with based on your comment @termi-official. Just to show what could be done, here are all three approaches. The third is the most involved but that's just how it is.

using OrderedCollections
function to_facetindex(cell_map, ijk, e)
    uvw = DT.sort_triangle(ijk)
    cell_number = cell_map[uvw]
    u, v, w = triangle_vertices(uvw)
    i, j = edge_vertices(e)
    (i, j) == (u, v) && return FacetIndex(cell_number, 1)
    (i, j) == (v, w) && return FacetIndex(cell_number, 2)
    return FacetIndex(cell_number, 3) # assume (i, j) == (w, u)
end
function get_facetsets_1(tri::Triangulation, cell_map)
    # Using each_boundary_edge (same as keys(get_boundary_edge_map(tri)))
    facetsets = Dict{String,OrderedSet{FacetIndex}}()
    ghost_map = Dict(i => string(i) for i in each_ghost_vertex(tri)) # trying to avoid allocating string(i) every single time
    for e in DT.each_boundary_edge(tri)
        u, v = edge_vertices(e)
        w = get_adjacent(tri, e)
        g = get_adjacent(tri, v, u) # ghost vertex
        set = get!(OrderedSet{FacetIndex}, facetsets, ghost_map[g])
        push!(set, to_facetindex(cell_map, (u, v, w), e))
    end
    return facetsets
end
function get_facetsets_2(tri::Triangulation, cell_map)
    # Using each_ghost_vertex together with get_adjacent2vertex 
    facetsets = Dict{String,OrderedSet{FacetIndex}}()
    for g in each_ghost_vertex(tri)
        edges = get_adjacent2vertex(tri, g)
        set = OrderedSet{FacetIndex}()
        facetsets[string(g)] = set
        for e′ in each_edge(edges)
            e = DT.reverse_edge(e′) # e′ is the edge adjoining the ghost vertex 
            u, v = edge_vertices(e)
            w = get_adjacent(tri, e)
            push!(set, to_facetindex(cell_map, (u, v, w), e))
        end
    end
    return facetsets
end
function get_facetsets_3(tri::Triangulation, cell_map)
    # Iterating over the boundary nodes in order
    facetsets = Dict{String,OrderedSet{FacetIndex}}()
    g = 0
    nc = DT.num_curves(tri)
    for curve_index in 1:nc
        if nc == 1 
            bn = get_boundary_nodes(tri)
        else
            bn = get_boundary_nodes(tri, curve_index)
        end
        ns = DT.num_sections(bn)
        for segment_index in 1:ns
            g -= 1
            set = OrderedSet{FacetIndex}()
            facetsets[string(g)] = set
            if nc == ns == 1
                bnn = bn 
            else
                bnn = get_boundary_nodes(bn, segment_index) 
            end
            ne = num_boundary_edges(bnn)
            for i in 1:ne
                u = get_boundary_nodes(bnn, i)
                v = get_boundary_nodes(bnn, i + 1)
                w = get_adjacent(tri, u, v)
                push!(set, to_facetindex(cell_map, (u, v, w), (u, v)))
            end
        end
    end
    return facetsets
end

A function to then generate a grid would be

function Ferrite.Grid(tri::Triangulation)
    cells, cell_map = get_cells(tri)
    nodes = get_nodes(tri)
    facetsets = get_facetsets_3(tri, cell_map)
    return Ferrite.Grid(cells, nodes; facetsets)
end

As an example, here I create a Grid for an annulus.

function Annulus(R₁, R₂) # centred at the origin 
    @assert R₁ < R₂
    inner = CircularArc((R₁, 0.0), (R₁, 0.0), (0.0, 0.0), positive = false)
    outer = CircularArc((R₂, 0.0), (R₂, 0.0), (0.0, 0.0))
    return [[[outer]], [[inner]]]
end
R₁, R₂ = 1.0, 2.0
tri = triangulate(NTuple{2, Float64}[], boundary_nodes = Annulus(R₁, R₂))
refine!(tri; max_area = 1e-3π * (R₂^2 - R₁^2))

julia> grid = Grid(tri)
Grid{2, Triangle, Float64} with 1614 Triangle cells and 890 nodes

julia> grid.facetsets
Dict{String, OrderedSet{FacetIndex}} with 2 entries:
  "-2" => OrderedSet{FacetIndex}([FacetIndex((1473, 3)), FacetIndex((129, 2)), FacetIndex((257, 3)), FacetIndex((196, 2)), FacetIndex((1534, 3)), FacetIndex((1295, 3)), FacetIndex((671, 2)), FacetInde…
  "-1" => OrderedSet{FacetIndex}([FacetIndex((686, 3)), FacetIndex((1469, 3)), FacetIndex((105, 2)), FacetIndex((331, 3)), FacetIndex((673, 2)), FacetIndex((252, 3)), FacetIndex((1463, 2)), FacetIndex…

What do you think?

I handle the boundary conditions a bit differently inside FiniteVolumeMethod.jl (https://github.com/SciML/FiniteVolumeMethod.jl/blob/9dc2612219b95eb49b4a81d8cace441f10b6ded4/src/conditions.jl#L452C1-L490C1), which also allows for internal constraints (useful for e.g. this example https://sciml.github.io/FiniteVolumeMethod.jl/stable/wyos/poissons_equation/#Using-the-Provided-Template). I don't know if you allow for internal constraints in this package-I imagine you do-, but it wouldn't be too hard for me to include that. Probably out of scope for this issue but just as an FYI I guess.

[DanielVandH]

Might be something wrong in my code above, not really sure. I tried using it for https://ferrite-fem.github.io/Ferrite.jl/dev/tutorials/heat_equation/, using triangles instead of quadrilaterals but I get some singular matrix:

using Ferrite, SparseArrays, DelaunayTriangulation
tri = triangulate_rectangle(-1, 1, -1, 1, 20, 20, single_boundary=true)
grid = Grid(tri)

ip = Lagrange{RefTriangle,1}()
qr = QuadratureRule{RefTriangle}(2)
cellvalues = CellValues(qr, ip)

dh = DofHandler(grid)
add!(dh, :u, ip)
close!(dh)

K = allocate_matrix(dh)

ch = ConstraintHandler(dh)
∂Ω = getfacetset(grid, "-1")

dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0.0)
add!(ch, dbc)
close!(ch)

function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues)
    n_basefuncs = getnbasefunctions(cellvalues)
    # Reset to 0
    fill!(Ke, 0)
    fill!(fe, 0)
    # Loop over quadrature points
    for q_point in 1:getnquadpoints(cellvalues)
        # Get the quadrature weight
        dΩ = getdetJdV(cellvalues, q_point)
        # Loop over test shape functions
        for i in 1:n_basefuncs
            δu = shape_value(cellvalues, q_point, i)
            ∇δu = shape_gradient(cellvalues, q_point, i)
            # Add contribution to fe
            fe[i] += δu * dΩ
            # Loop over trial shape functions
            for j in 1:n_basefuncs
                ∇u = shape_gradient(cellvalues, q_point, j)
                # Add contribution to Ke
                Ke[i, j] += (∇δu ⋅ ∇u) * dΩ
            end
        end
    end
    return Ke, fe
end

function assemble_global(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler)
    # Allocate the element stiffness matrix and element force vector
    n_basefuncs = getnbasefunctions(cellvalues)
    Ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)
    # Allocate global force vector f
    f = zeros(ndofs(dh))
    # Create an assembler
    assembler = start_assemble(K, f)
    # Loop over all cels
    for cell in CellIterator(dh)
        # Reinitialize cellvalues for this cell
        reinit!(cellvalues, cell)
        # Compute element contribution
        assemble_element!(Ke, fe, cellvalues)
        # Assemble Ke and fe into K and f
        assemble!(assembler, celldofs(cell), Ke, fe)
    end
    return K, f
end

K, f = assemble_global(cellvalues, K, dh)

ERROR: ArgumentError: det(J) is not positive: det(J) = -1.35180055401662
Stacktrace:
 [1] throw_detJ_not_pos(detJ::Float64)
   @ Ferrite C:\Users\danjv\.julia\packages\Ferrite\Hp5x1\src\FEValues\common_values.jl:5
 [2] _update_detJdV!
   @ C:\Users\danjv\.julia\packages\Ferrite\Hp5x1\src\FEValues\CellValues.jl:96 [inlined]
 [3] reinit!(cv::CellValues{…}, cell::Nothing, x::Vector{…})
   @ Ferrite C:\Users\danjv\.julia\packages\Ferrite\Hp5x1\src\FEValues\CellValues.jl:119
 [4] reinit!
   @ C:\Users\danjv\.julia\packages\Ferrite\Hp5x1\src\FEValues\CellValues.jl:102 [inlined]
 [5] reinit!
   @ C:\Users\danjv\.julia\packages\Ferrite\Hp5x1\src\iterators.jl:92 [inlined]
 [6] assemble_global(cellvalues::CellValues{…}, K::SparseMatrixCSC{…}, dh::DofHandler{…})
   @ Main .\Untitled-1:186
 [7] top-level scope
   @ Untitled-1:195
Some type information was truncated. Use `show(err)` to see complete types.

Will try look some more into it, but just posting for now in case there's anything obvious I'm missing.

[fredrikekre]

That is probably because Ferrite requires the node order to be counter clockwise. (Just skimming through, I guess the sorting of the node ids (DT.sort_triangl) could be the problem?)

[DanielVandH]

sort_triangle preserves the order of the triangle (counter-clockwise), it just rotates it to put the minimum vertex first. I think the orientation is consistent throughout but I'll check.

for T in grid.cells
    ijk = T.nodes
    DT.is_positively_oriented(DT.triangle_orientation(tri, ijk)) || error("Triangle is negatively oriented")
end

shows it's fine.

[DanielVandH]

Nevermind. It's because I fixed an issue in my post above but somehow didn't do the same in my actual script 🤦 .

[DanielVandH]

Here's a complete example. It's https://ferrite-fem.github.io/Ferrite.jl/dev/tutorials/heat_equation/ except on an annulus. The boundary conditions are zero on the inner boundary, and $x^2y^2$ on the outer boundary. I use tricontourf to plot the solution.

using Ferrite, SparseArrays, CairoMakie, DelaunayTriangulation, OrderedCollections

const DT = DelaunayTriangulation
function get_cells(tri::Triangulation)
    cells = Vector{Triangle}(undef, num_solid_triangles(tri))
    cell_map = Dict{NTuple{3,Int},Int}()
    for (cell_number, triangle) in enumerate(each_solid_triangle(tri))
        ijk = triangle_vertices(triangle)
        uvw = DT.sort_triangle(ijk) # It's useful to know the order of the vertices
        cells[cell_number] = Triangle(uvw)
        cell_map[uvw] = cell_number
    end
    return cells, cell_map
end
function get_nodes(tri::Triangulation)
    nodes = Vector{Node{2,Float64}}(undef, DT.num_points(tri))
    for vertex in DT.each_point_index(tri)
        p = get_point(tri, vertex)
        x, y = getxy(p)
        nodes[vertex] = Node((x, y))
    end
    return nodes
end

function to_facetindex(cell_map, ijk, e)
    uvw = DT.sort_triangle(ijk)
    cell_number = cell_map[uvw]
    u, v, w = triangle_vertices(uvw)
    i, j = edge_vertices(e)
    (i, j) == (u, v) && return FacetIndex(cell_number, 1)
    (i, j) == (v, w) && return FacetIndex(cell_number, 2)
    return FacetIndex(cell_number, 3) # assume (i, j) == (w, u)
end
function get_facetsets(tri::Triangulation, cell_map)
    # Iterating over the boundary nodes in order
    facetsets = Dict{String,OrderedSet{FacetIndex}}()
    g = 0
    nc = DT.num_curves(tri)
    for curve_index in 1:nc
        if nc == 1
            bn = get_boundary_nodes(tri)
        else
            bn = get_boundary_nodes(tri, curve_index)
        end
        ns = DT.num_sections(bn)
        for segment_index in 1:ns
            g -= 1
            set = OrderedSet{FacetIndex}()
            facetsets[string(g)] = set
            if nc == ns == 1
                bnn = bn
            else
                bnn = get_boundary_nodes(bn, segment_index)
            end
            ne = num_boundary_edges(bnn)
            for i in 1:ne
                u = get_boundary_nodes(bnn, i)
                v = get_boundary_nodes(bnn, i + 1)
                w = get_adjacent(tri, u, v)
                push!(set, to_facetindex(cell_map, (u, v, w), (u, v)))
            end
        end
    end
    return facetsets
end

function Ferrite.Grid(tri::Triangulation)
    cells, cell_map = get_cells(tri)
    nodes = get_nodes(tri)
    facetsets = get_facetsets(tri, cell_map)
    return Ferrite.Grid(cells, nodes; facetsets)
end

function Annulus(R₁, R₂) # centred at the origin 
    @assert R₁ < R₂
    inner = CircularArc((R₁, 0.0), (R₁, 0.0), (0.0, 0.0), positive=false)
    outer = CircularArc((R₂, 0.0), (R₂, 0.0), (0.0, 0.0))
    return [[[outer]], [[inner]]]
end
R₁, R₂ = 1.0, 2.0
tri = triangulate(NTuple{2,Float64}[], boundary_nodes=Annulus(R₁, R₂))
refine!(tri; max_area=1e-3π * (R₂^2 - R₁^2))
grid = Grid(tri)

ip = Lagrange{RefTriangle,1}()
qr = QuadratureRule{RefTriangle}(2)
cellvalues = CellValues(qr, ip)
dh = DofHandler(grid)
add!(dh, :u, ip)
close!(dh)
K = allocate_matrix(dh)

ch = ConstraintHandler(dh)
dbc = Dirichlet(:u, getfacetset(grid, "-1"), (x, t) -> x[1]^2 * x[2]^2)
add!(ch, dbc)
nbc = Dirichlet(:u, getfacetset(grid, "-2"), (x, t) -> 0.0)
add!(ch, nbc)
close!(ch)

function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellValues)
    n_basefuncs = getnbasefunctions(cellvalues)
    # Reset to 0
    fill!(Ke, 0)
    fill!(fe, 0)
    # Loop over quadrature points
    for q_point in 1:getnquadpoints(cellvalues)
        # Get the quadrature weight
        dΩ = getdetJdV(cellvalues, q_point)
        # Loop over test shape functions
        for i in 1:n_basefuncs
            δu = shape_value(cellvalues, q_point, i)
            ∇δu = shape_gradient(cellvalues, q_point, i)
            # Add contribution to fe
            fe[i] += δu * dΩ
            # Loop over trial shape functions
            for j in 1:n_basefuncs
                ∇u = shape_gradient(cellvalues, q_point, j)
                # Add contribution to Ke
                Ke[i, j] += (∇δu ⋅ ∇u) * dΩ
            end
        end
    end
    return Ke, fe
end

function assemble_global(cellvalues::CellValues, K::SparseMatrixCSC, dh::DofHandler)
    # Allocate the element stiffness matrix and element force vector
    n_basefuncs = getnbasefunctions(cellvalues)
    Ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)
    # Allocate global force vector f
    f = zeros(ndofs(dh))
    # Create an assembler
    assembler = start_assemble(K, f)
    # Loop over all cels
    for cell in CellIterator(dh)
        # Reinitialize cellvalues for this cell
        reinit!(cellvalues, cell)
        # Compute element contribution
        assemble_element!(Ke, fe, cellvalues)
        # Assemble Ke and fe into K and f
        assemble!(assembler, celldofs(cell), Ke, fe)
    end
    return K, f
end

K, f = assemble_global(cellvalues, K, dh)
apply!(K, f, ch)
u = K \ f
points = [p.x for p in getnodes(grid)]
projector = L2Projector(ip, grid)
ph = PointEvalHandler(grid, points)
sol = evaluate_at_points(ph, projector, u)
tricontourf(tri, sol)

Let me know if there is anything else you think could be useful to demonstrate.

[KnutAM]

That looks very nice! If you want to, it would be great to open a PR adding this as a how-to IMO. Then it is also easy to comment / test it!

Regarding the postprocessing, I'd suggest

- points = [p.x for p in getnodes(grid)]
- projector = L2Projector(ip, grid)
- ph = PointEvalHandler(grid, points)
- sol = evaluate_at_points(ph, projector, u)
+ sol =  evaluate_at_grid_nodes(dh, u, :u)

(The fact that evaluate_at_points(ph, projector, u) works is a bit unfortunate in this case actually, since u then should be the result of a projection and be numbered according to the projector's DofHandler, but I don't know how we could detect this. In this case, it will give the correct results since the dof handlers are identical.)

[DanielVandH]

Thanks for the tip with evaluate_at_grid_nodes! I figured there was an easier way but that was the best I could find in the docs.

I'll write up the example properly with Literate and make a PR soon

[DanielVandH]

See #1040