README.md

PR 780 P1 reproductions

These standalone scripts reproduce the three P1 findings from the review. Run them from the Ferrite.jl repository root. Each script currently fails at the documented invariant; after a fix, it should complete without error.

julia --project=. review_mwes/p1_balanceforest_3d_corner.jl
julia --project=. review_mwes/p1_q2_conformity.jl
julia --project=. review_mwes/p1_subdomain_conformity.jl

• p1_balanceforest_3d_corner.jl: balanceforest! leaves two corner-touching 3D leaves at levels 3 and 1.
• p1_q2_conformity.jl: the Q1 averaging constraint is silently applied to a Q2 field, so an exactly representable quadratic is not reproduced.
• p1_subdomain_conformity.jl: a projector restricted to a subdomain creates a hanging constraint for absent dof zero and fails while closing its constraint handler.

p1_balanceforest_3d_corner.jl

using Ferrite

function leaf_at_corner(tree, corner)
    return only(o for o in tree.leaves if corner in Ferrite.AMR.vertices(o, tree.b))
end

forest = ForestBWG(generate_grid(Hexahedron, (2, 2, 2)), 7)

# Trees 4 and 8 meet at one macro-grid corner. Refine tree 4 three times toward
# that corner, then ask the forest to restore its documented full 2:1 balance.
tree4 = forest.cells[4]
corner4 = (0, 0, 1 << tree4.b)
for _ in 1:3
    Ferrite.refine!(tree4, leaf_at_corner(tree4, corner4))
end
Ferrite.balanceforest!(forest)

tree8 = forest.cells[8]
corner8 = (1 << tree8.b, 1 << tree8.b, 0)
level4 = Int(leaf_at_corner(tree4, corner4).l)
level8 = Int(leaf_at_corner(tree8, corner8).l)

println("levels at the shared corner: tree 4 = $level4, tree 8 = $level8")
@assert abs(level4 - level8) <= 1 "balanceforest! left a 3:1 corner connection"

p1_q2_conformity.jl

using Ferrite

function reproduce_q2_projection_error()
    forest = ForestBWG(generate_grid(Quadrilateral, (2, 1)), 3)
    Ferrite.refine!(forest, [1])
    Ferrite.balanceforest!(forest)
    grid = Ferrite.creategrid(forest)

    ip = Lagrange{RefQuadrilateral, 2}()
    qr = QuadratureRule{RefQuadrilateral}(4)
    cv = CellValues(qr, ip, geometric_interpolation(Quadrilateral))

    projector = L2Projector(grid)
    add!(projector, collect(1:getncells(grid)), ip; qr_rhs = qr)
    close!(projector)

    # y^2 is exactly representable by Q2. The projection should reproduce it to
    # roundoff, including across the hanging interface.
    qp_data = map(1:getncells(grid)) do cellid
        coordinates = getcoordinates(grid, cellid)
        reinit!(cv, coordinates)
        map(1:getnquadpoints(cv)) do qp
            y = spatial_coordinate(cv, qp, coordinates)[2]
            y^2
        end
    end
    projected = project(projector, qp_data, qr)

    max_error = 0.0
    for cell in CellIterator(projector.dh)
        coordinates = getcoordinates(cell)
        reinit!(cv, coordinates)
        element_values = projected[celldofs(cell)]
        for qp in 1:getnquadpoints(cv)
            y = spatial_coordinate(cv, qp, coordinates)[2]
            max_error = max(max_error, abs(function_value(cv, qp, element_values) - y^2))
        end
    end
    return max_error
end

max_error = reproduce_q2_projection_error()
println("maximum Q2 projection error: $max_error")
@assert max_error <= 1.0e-10 "ConformityConstraint applied Q1 averaging weights to Q2"

p1_subdomain_conformity.jl

using Ferrite

forest = ForestBWG(generate_grid(Quadrilateral, (2, 1)), 3)
Ferrite.refine!(forest, [1])
Ferrite.balanceforest!(forest)
grid = Ferrite.creategrid(forest)

# Refining the first macro cell produces cells 1:4; cell 5 is the unrefined
# coarse macro cell. The projector API supports selecting only that subdomain.
@assert getncells(grid) == 5
projector = L2Projector(grid)
add!(projector, Set([5]), Lagrange{RefQuadrilateral, 1}(); qr_rhs = nothing)

# Current result: BoundsError indexing a BitVector at dof zero while closing the
# automatically created conformity ConstraintHandler.
close!(projector)
println(projector)