Plotting surface for second order triangular mesh fails

[mipals]

Hey there,

As discussed on slack, here is a MWE that shows that surface fails to plot directly when the geometry is quadratic with the following error

ERROR: ArgumentError: the number of (geometric) base functions (3) does not match the number of coordinates in the vector (6). Perhaps you forgot to use an appropriate geometric interpolation when creating FE values? See https://github.com/Ferrite-FEM/Ferrite.jl/issues/265 for more details.
Stacktrace:
  [1] throw_incompatible_coord_length(length_x::Int64, n_base_funcs::Int64)
    @ Ferrite ~/.julia/packages/Ferrite/vFoCZ/src/FEValues/common_values.jl:20
  [2] reinit!
    @ ~/.julia/packages/Ferrite/vFoCZ/src/FEValues/CellValues.jl:115 [inlined]
  [3] reinit!
    @ ~/.julia/packages/Ferrite/vFoCZ/src/FEValues/CellValues.jl:102 [inlined]
  [4] reinit!(pv::PointValues{CellValues{…}}, x::Vector{Tensors.Vec{…}}, ξ::Tensors.Vec{2, Float64})
    @ Ferrite ~/.julia/packages/Ferrite/vFoCZ/src/FEValues/PointValues.jl:74
  [5] _transfer_solution!(data::Matrix{…}, pv::PointValues{…}, sdh::SubDofHandler{…}, ip_field::Lagrange{…}, cellset_::Vector{…}, val_buffer::Float64, val::Float64, field_name::Symbol, field_dim::Int64, plotter::MakiePlotter{…}, u::Vector{…}, process::typeof(FerriteViz.postprocess))
    @ FerriteViz ~/.julia/packages/FerriteViz/bXnNg/src/utils.jl:412
  [6] transfer_solution(plotter::MakiePlotter{…}, u::Vector{…}; field_name::Symbol, process::typeof(FerriteViz.postprocess))
    @ FerriteViz ~/.julia/packages/FerriteViz/bXnNg/src/utils.jl:389
  [7] (::FerriteViz.var"#103#105"{MakiePlotter{…}})(arg1#270::Vector{Float64}, arg2#271::Symbol, arg3#272::Function)
    @ FerriteViz ~/.julia/packages/FerriteViz/bXnNg/src/makieplotting.jl:327
  [8] #map#13
    @ ~/.julia/packages/Observables/YdEbO/src/Observables.jl:570 [inlined]
  [9] map(::FerriteViz.var"#103#105"{MakiePlotter{…}}, ::Observable{Vector{…}}, ::Observable{Any}, ::Observable{Any})
    @ Observables ~/.julia/packages/Observables/YdEbO/src/Observables.jl:568
 [10] plot!(SF::Plot{FerriteViz.surface, Tuple{MakiePlotter{…}}})
    @ FerriteViz ~/.julia/packages/FerriteViz/bXnNg/src/makieplotting.jl:324
 [11] connect_plot!(parent::Scene, plot::Plot{FerriteViz.surface, Tuple{MakiePlotter{…}}})
    @ Makie ~/.julia/packages/Makie/eERNK/src/interfaces.jl:395
 [12] plot!
    @ ~/.julia/packages/Makie/eERNK/src/interfaces.jl:412 [inlined]
 [13] plot!(ax::LScene, plot::Plot{FerriteViz.surface, Tuple{MakiePlotter{…}}})
    @ Makie ~/.julia/packages/Makie/eERNK/src/figureplotting.jl:412
 [14] plot!(fa::Makie.FigureAxis, plot::Plot{FerriteViz.surface, Tuple{MakiePlotter{…}}})
    @ Makie ~/.julia/packages/Makie/eERNK/src/figureplotting.jl:407
 [15] _create_plot(F::Function, attributes::Dict{…}, args::MakiePlotter{…})
    @ Makie ~/.julia/packages/Makie/eERNK/src/figureplotting.jl:318
 [16] surface(args::MakiePlotter{…})
    @ FerriteViz ~/.julia/packages/MakieCore/rTINf/src/recipes.jl:188

First of all the error is a little misleading, as it is not the user arguments that is wrong, but rather something internally in surface.

I can see the benefit of not directly plotting second order geometry as linear approximations, as in theory one would be plotting the "wrong" thing. However, wireframe does exactly this so supporting the same for surface seems reasonable.

In the MWE i did include how I ended up plotting what I wanted by extracting the vertices of the element and then creating a new linear grid and corresponding DofHandler.

MWE

using FerriteGmsh
using FerriteViz #do/ferrite-1.0-prep
using Ferrite
using KrylovKit
using GLMakie

function create_disk(refinements;order=1,a=1.0,R=3.0)
    @assert a < R && (a > 0) && (R > 0) && (order < 3) && (order > 0)
    gmsh.initialize()
    gmsh.model.occ.addDisk(0.0,0.0,0.0,a,a,1)
    gmsh.model.occ.addDisk(0.0,0.0,0.0,R,R,2)
    gmsh.model.occ.cut([(2,2)],[(2,1)])
    gmsh.model.occ.synchronize()
    gmsh.model.mesh.generate(2)
    # To get good solution quality refine the elements several times
    for _ in 1:refinements
        gmsh.model.mesh.refine()
    end
    gmsh.model.mesh.setOrder(order)
    nodes = tonodes()
    elements, _ = toelements(2)
    gmsh.finalize()
    grid = Grid(elements, nodes);
    return grid
end

# Creating the disk
refinements = 0
order = 2
grid = create_disk(refinements;order=order)

FerriteViz.wireframe(grid)

# Defining the FEM assembly
function assemble_elements!(Ke, Me, cellvalues::CellValues)
    n_basefuncs = getnbasefunctions(cellvalues)
    ## Reset to 0
    fill!(Ke, 0)
    fill!(Me, 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)
            ## Loop over trial shape functions
            for j in 1:n_basefuncs
                u = shape_value(cellvalues, q_point, j)
                ∇u = shape_gradient(cellvalues, q_point, j)
                ## Adding contributions to element matrices
                Ke[i, j] += (∇δu ⋅ ∇u) * dΩ
                Me[i, j] += (δu ⋅ u) * dΩ
            end
        end
    end
    return Ke, Me
end;

# Now what is left is to define a function that sums up all elements contributions
function assemble_global(cellvalues::CellValues, K, M, dh::DofHandler)
    ## Allocate the element stiffness and mass matrices
    n_basefuncs = getnbasefunctions(cellvalues)
    Me = zeros(n_basefuncs, n_basefuncs)
    Ke = zeros(n_basefuncs, n_basefuncs)
    ## Create assemblers for both M and K
    assemblerM = start_assemble(M) 
    assemblerK = start_assemble(K) 
    ## Loop over all cels
    for cell in CellIterator(dh)
        ## Reinitialize cellvalues for this cell
        reinit!(cellvalues, cell)
        ## Compute element contribution
        assemble_elements!(Ke, Me, cellvalues)
        ## Insert Ke and Me into their respective matrices
        assemble!(assemblerM, celldofs(cell), Me)
        assemble!(assemblerK, celldofs(cell), Ke)
    end
    return K, M
end;

# Defining the interpolation
ip = Lagrange{RefTriangle, 2}()
ip_geo = Lagrange{RefTriangle, order}()
qr = QuadratureRule{RefTriangle}(2)
cellvalues = CellValues(qr, ip, ip_geo)

# Creating the DofHandler and getting the pressure
dh = DofHandler(grid)
add!(dh, :p, ip)
close!(dh)

# Allocating stiffness and mass matrix
K = allocate_matrix(dh)
M = allocate_matrix(dh);

# Assembling stiffness and mass matrices
@time K,M = assemble_global(cellvalues, K, M, dh);

# In the case of a triangle the there exist an analytical solution to the problem
λ, X, conv = geneigsolve((K,M), dh.ndofs, 10, :SR,  issymmetric=true, isposdef=true,maxiter=400)

# Trying to plot the first non-zero eigenvector as
plotter = FerriteViz.MakiePlotter(dh,X[2])
FerriteViz.wireframe(plotter) # Works
FerriteViz.surface(plotter)   # Does not work

# Casting to a linear grid seem to resolve the issue?
lin_cells = [Triangle(cell.nodes[1:3]) for cell in grid.cells]
lin_grid = Grid(lin_cells, grid.nodes) # Keeping unused nodes
lin_dh = DofHandler(lin_grid)
add!(lin_dh, :p, ip)
close!(lin_dh)
lin_plotter = FerriteViz.MakiePlotter(lin_dh,X[2])
FerriteViz.wireframe(lin_plotter) # Works. Unused nodes still shown?
FerriteViz.surface(lin_plotter)

[fredrikekre]

I pushed a commit that fixes so try again after pkg> update FerriteViz (see 041b32c)

[mipals]

That does seem to work! Thanks 👍

[koehlerson]

Thanks @fredrikekre

Fixed in #103