[X-PR] p4est edge operations and consistent corner operation

[koehlerson]

• intraoctree
  • edge_neighbor
  • tests
  • docs
• interoctree
  • remote version that implements exactly the logic from the paper
  • version with flipped logic that seems more natural (analogous to transform_face)
  • tests
  • docs
• high level functionality
  • balancing
  • creategrid
  • hanging nodes

[codecov]

Codecov Report

Attention: Patch coverage is 94.76190% with 11 lines in your changes are missing coverage. Please review.

❗ No coverage uploaded for pull request base (mk/p4est@0611729). Click here to learn what that means.

❗ Current head 160d0d9 differs from pull request most recent head f69cc60. Consider uploading reports for the commit f69cc60 to get more accurate results

Files	Patch %	Lines
src/Adaptivity/AdaptiveCells.jl	95.21%	10 Missing ⚠️
src/Export/VTK.jl	0.00%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##             mk/p4est     #902   +/-   ##
===========================================
  Coverage            ?   33.50%           
===========================================
  Files               ?       36           
  Lines               ?     6049           
  Branches            ?        0           
===========================================
  Hits                ?     2027           
  Misses              ?     4022           
  Partials            ?        0           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[termi-official]

Suggested change

	function OctantBWG(level::T,coords::NTuple) where T <: Integer
	function OctantBWG(level::Integer,coords::NTuple)

Is there a reason why you grab T explicitly?

[koehlerson]

I thought that we could add later a check if level and eltype(coords) matches and use the wider/smaller type

[termi-official]

I think we should start adding some docs to the code to understand what we have created here :D

[termi-official]

Seems like rebalance over edges is still broken. MWE:

using Ferrite, FerriteGmsh, SparseArrays
grid = generate_grid(Hexahedron, (3,3,2));
grid  = ForestBWG(grid,10)

function adapt_mesh_random_sequence(initial_grid, nctor = 5)
    i = 1
    grid = deepcopy(initial_grid)
    pvd = paraview_collection("amr-grid.pvd");
    while i<=10
        @show i
        transfered_grid = Ferrite.creategrid(grid)
        vtk_grid("amr-grid_$i", transfered_grid) do vtk
            pvd[i] = vtk
        end

        dh = DofHandler(transfered_grid)
        add!(dh, :u, Lagrange{RefHexahedron, 1}())
        close!(dh);

        ch = ConstraintHandler(dh)
        add!(ch, ConformityConstraint(:u))
        close!(ch);

        cells_to_refine = rand(1:length(transfered_grid.cells), nctor)
        Ferrite.refine!(grid, cells_to_refine)
        Ferrite.balanceforest!(grid)

        i += 1
    end
    vtk_save(pvd);
end

adapt_mesh_random_sequence(grid)

Edit: Deterministic reproducer:

using Ferrite, FerriteGmsh, SparseArrays
grid = generate_grid(Hexahedron, (2,1,1));
grid  = ForestBWG(grid,10)

function adapt_mesh_random_sequence(initial_grid, nctor = 5)
    grid = deepcopy(initial_grid)

    transfered_grid = Ferrite.creategrid(grid)
    vtk_grid("amr-grid_initial", transfered_grid) do vtk
    end

    Ferrite.refine!(grid, [1,2])
    Ferrite.balanceforest!(grid)
    Ferrite.refine!(grid, [4])
    Ferrite.balanceforest!(grid)
    Ferrite.refine!(grid, [5])
    Ferrite.balanceforest!(grid)
    
    transfered_grid = Ferrite.creategrid(grid)
    vtk_grid("amr-grid_broken", transfered_grid) do vtk
    end
end

adapt_mesh_random_sequence(grid)

[koehlerson]

looks correctly balanced to me, but throws nested affine constraints

[termi-official]

Dumping this one for later

using Ferrite, FerriteGmsh, SparseArrays
grid = generate_grid(Hexahedron, (3,3,2));
function random_deformation_field(x::Vec{sdim}) where sdim
    if any(x .≈ -1.0) || any(x .≈ 1.0)
        return x
    else
        Vec{sdim}(x .+ (rand(sdim).-0.5)*0.15)
    end
end
transform_coordinates!(grid, random_deformation_field)
grid  = ForestBWG(grid,10)

analytical_solution(x) = atan(2*(norm(x)-0.5)/0.02)
analytical_rhs(x) = -laplace(analytical_solution,x)

function assemble_cell!(ke, fe, cellvalues, ue, coords)
    fill!(ke, 0.0)
    fill!(fe, 0.0)

    n_basefuncs = getnbasefunctions(cellvalues)
    for q_point in 1:getnquadpoints(cellvalues)
        x = spatial_coordinate(cellvalues, q_point, coords)
        dΩ = getdetJdV(cellvalues, q_point)
        for i in 1:n_basefuncs
            Nᵢ = shape_value(cellvalues, q_point, i)
            ∇Nᵢ = shape_gradient(cellvalues, q_point, i)# shape_symmetric_gradient(cellvalues, q_point, i)
            fe[i] += analytical_rhs(x) * Nᵢ * dΩ # add internal force to residual
            for j in 1:n_basefuncs
                ∇Nⱼ = shape_gradient(cellvalues, q_point, j)
                ke[i, j] += ∇Nⱼ ⋅ ∇Nᵢ * dΩ
            end
        end
    end
end

function assemble_global!(K, f, a, dh, cellvalues)
    ## Allocate the element stiffness matrix and element force vector
    n_basefuncs = getnbasefunctions(cellvalues)
    ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)
    ## Create an assembler
    assembler = start_assemble(K, f)
    ## Loop over all cells
    for cell in CellIterator(dh)
        reinit!(cellvalues, cell) # update spatial derivatives based on element coordinates
        @views ue = a[celldofs(cell)]
        ## Compute element contribution
        coords = getcoordinates(cell)
        assemble_cell!(ke, fe, cellvalues, ue, coords)
        ## Assemble ke and fe into K and f
        assemble!(assembler, celldofs(cell), ke, fe)
    end
    return K, f
end

function solve(grid)
    dim = 2
    order = 1 # linear interpolation
    ip = Lagrange{RefHexahedron, order}() # vector valued interpolation
    qr = QuadratureRule{RefHexahedron}(2) # 1 quadrature point
    cellvalues = CellValues(qr, ip);

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

    ch = ConstraintHandler(dh)
    add!(ch, ConformityConstraint(:u))
    add!(ch, Dirichlet(:u, getfaceset(grid, "top"), (x, t) -> 0.0))
    add!(ch, Dirichlet(:u, getfaceset(grid, "right"), (x, t) -> 0.0))
    add!(ch, Dirichlet(:u, getfaceset(grid, "left"), (x, t) -> 0.0))
    add!(ch, Dirichlet(:u, getfaceset(grid, "bottom"), (x, t) -> 0.0))
    close!(ch);

    K = create_sparsity_pattern(dh,ch)
    f = zeros(ndofs(dh))
    a = zeros(ndofs(dh))
    assemble_global!(K, f, a, dh, cellvalues);
    apply!(K, f, ch)
    u = K \ f;
    apply!(u,ch)
    return u,dh,ch,cellvalues
end

function compute_fluxes(u,dh::DofHandler{sdim}) where sdim
    ip = Lagrange{RefHexahedron, 1}()
    # Normal quadrature points
    qr = QuadratureRule{RefHexahedron}(2)
    cellvalues = CellValues(qr, ip);
    # Superconvergent point
    qr_sc = QuadratureRule{RefHexahedron}(1)
    cellvalues_sc = CellValues(qr_sc, ip);
    #Buffers
    σ_gp = Vector{Vector{Vec{sdim,Float64}}}()
    σ_gp_loc = Vector{Vec{sdim,Float64}}()
    σ_gp_sc = Vector{Vector{Vec{sdim,Float64}}}()
    σ_gp_sc_loc = Vector{Vec{sdim,Float64}}()
    for (cellid,cell) in enumerate(CellIterator(dh))
        @views ue = u[celldofs(cell)]

        reinit!(cellvalues, cell)
        for q_point in 1:getnquadpoints(cellvalues)
            gradu = function_gradient(cellvalues, q_point, ue)
            push!(σ_gp_loc, gradu)
        end
        push!(σ_gp,copy(σ_gp_loc))
        empty!(σ_gp_loc)

        reinit!(cellvalues_sc, cell)
        for q_point in 1:getnquadpoints(cellvalues_sc)
            gradu = function_gradient(cellvalues_sc, q_point, ue)
            push!(σ_gp_sc_loc, gradu)
        end
        push!(σ_gp_sc,copy(σ_gp_sc_loc))
        empty!(σ_gp_sc_loc)
    end
    return σ_gp, σ_gp_sc
end

function solve_adaptive(initial_grid::Ferrite.AbstractGrid{sdim}) where sdim
    ip = Lagrange{RefHexahedron, 1}()^sdim
    qr_sc = QuadratureRule{RefHexahedron}(1)
    cellvalues_flux = CellValues(qr_sc, ip);
    finished = false
    i = 1
    grid = deepcopy(initial_grid)
    pvd = paraview_collection("heat_amr.pvd");
    while !finished && i<=10
        @show i
        transfered_grid = Ferrite.creategrid(grid)
        vtk_grid("heat_amr-grid_$i", transfered_grid) do vtk
        end
        u,dh,ch,cv = solve(transfered_grid)
        σ_gp, σ_gp_sc = compute_fluxes(u,dh)
        projector = L2Projector(Lagrange{RefHexahedron, 1}(), transfered_grid)
        σ_dof = project(projector, σ_gp, QuadratureRule{RefHexahedron}(2))
        cells_to_refine = Int[]
        error_arr = Float64[]
        for (cellid,cell) in enumerate(CellIterator(projector.dh))
            reinit!(cellvalues_flux, cell)
            @views σe = σ_dof[celldofs(cell)]
            error = 0.0
            for q_point in 1:getnquadpoints(cellvalues_flux)
                σ_dof_at_sc = function_value(cellvalues_flux, q_point, σe)
                error += norm((σ_gp_sc[cellid][q_point] - σ_dof_at_sc ))
                error *= getdetJdV(cellvalues_flux,q_point)
            end
            if error > 0.001
                push!(cells_to_refine,cellid)
            end
            push!(error_arr,error)
        end

        vtk_grid("heat_amr-iteration_$i", dh) do vtk
            vtk_point_data(vtk, dh, u)
            vtk_point_data(vtk, projector, σ_dof, "flux")
            vtk_cell_data(vtk, getindex.(collect(Iterators.flatten(σ_gp_sc)),1), "flux sc x")
            vtk_cell_data(vtk, getindex.(collect(Iterators.flatten(σ_gp_sc)),2), "flux sc y")
            vtk_cell_data(vtk, error_arr, "error")
            pvd[i] = vtk
        end

        Ferrite.refine!(grid, cells_to_refine)
        Ferrite.balanceforest!(grid)

        i += 1
        if isempty(cells_to_refine)
            finished = true
        end
    end
    vtk_save(pvd);
end

solve_adaptive(grid)