begin
    using ExtendableFEMBase
    using ExtendableGrids
    using ExtendableSparse
    using GridVisualize
    using PlutoVista
    GridVisualize.default_plotter!(PlutoVista)
end

PlutoVista

begin
    ## PDE data
    μ = 1.0
    f = x -> x[1] - x[2]

    ## discretization parameters
    nref = 9
    order = 1
end

1

tricontour(xgrid[Coordinates],xgrid[CellNodes],sol.entries[1:num_nodes(xgrid)]; levels = 5, resolution = (500,500))

begin
    ## call low level solver
    sol = solve_poisson_lowlevel(FES, μ, f)
end

FEVector information
====================
   block  |  ndofs 	|     min  /  max    	| FEType 		 (name/tag)
 [    1]  |  263169	| -1.21e-02/1.21e-02  	| H1Pk{1,2,1}  	 (#1)

function solve_poisson_lowlevel(FES, μ, f)
    
    Solution = FEVector(FES)
    FES = Solution[1].FES
    A = FEMatrix(FES, FES)
    b = FEVector(FES)
    println("Assembling operators...")
    @time assemble!(A.entries, b.entries, FES, f, μ)

    ## fix boundary dofs
    println("Assembling boundary data...")
    @time begin
        BFaceDofs::Adjacency{Int32} = FES[ExtendableFEMBase.BFaceDofs]
        nbfaces::Int = num_sources(BFaceDofs)
        AM::ExtendableSparseMatrix{Float64,Int64} = A.entries
        dof_j::Int = 0
        for bface = 1 : nbfaces
            for j = 1 : num_targets(BFaceDofs,1)
                dof_j = BFaceDofs[j, bface]
                AM[dof_j,dof_j] = 1e60
                b.entries[dof_j] = 0
            end
        end
    end
    ExtendableSparse.flush!(A.entries)

    ## solve
    println("Solving linear system...")
    @time copyto!(Solution.entries, A.entries \ b.entries)

    return Solution
end

solve_poisson_lowlevel (generic function with 1 method)

begin
    ## create finite element space
    FEType = H1Pk{1,2,order}

    ## prepare finite element space and dofmaps
    println("Creating FESpace...")
    @time FES = FESpace{FEType}(xgrid)
    FES
end

FESpace information
===================
     name = H1Pk{1,2,1}
   FEType = H1Pk{1,2,1}
  FEClass = ExtendableFEMBase.AbstractH1FiniteElement
    ndofs = 263169


DofMaps
==========

begin
    ## create grid
    X = LinRange(0,1,2^nref+1)
    Y = LinRange(0,1,2^nref+1)
    println("Creating grid...")
    @time xgrid = simplexgrid(X,Y)
    println("Preparing FaceNodes...")
    @time xgrid[FaceNodes]
    println("Preparing CellVolumes...")
    @time xgrid[CellVolumes]
    xgrid
end

ExtendableGrids.ExtendableGrid{Float64, Int32}
      dim =       2
   nnodes =  263169
   ncells =  524288
  nbfaces =    2048

function assemble!(A::ExtendableSparseMatrix, b::Vector, FES, f, μ = 1)

    xgrid = FES.xgrid
    EG = xgrid[UniqueCellGeometries][1]
    FEType = eltype(FES)
    L2G = L2GTransformer(EG, xgrid, ON_CELLS)

    ## dofmap
    CellDofs = FES[ExtendableFEMBase.CellDofs]
    
    ## quadrature formula
    qf = QuadratureRule{Float64, EG}(2*(get_polynomialorder(FEType, EG)-1))
    weights::Vector{Float64} = qf.w
    xref::Vector{Vector{Float64}} = qf.xref
    nweights::Int = length(weights)
    
    ## FE basis evaluator
    FEBasis_∇ = FEEvaluator(FES, Gradient, qf)
    ∇vals = FEBasis_∇.cvals
    FEBasis_id = FEEvaluator(FES, Identity, qf)
    idvals = FEBasis_id.cvals

    
    cellvolumes = xgrid[CellVolumes]
   
    ## ASSEMBLY LOOP
    function barrier(EG, L2G::L2GTransformer)
        ## barrier function to avoid allocations by EG dispatch
        
    	ndofs4cell::Int = get_ndofs(ON_CELLS, FEType, EG)
    	Aloc = zeros(Float64, ndofs4cell, ndofs4cell)
        ncells::Int = num_cells(xgrid)
    	dof_j::Int, dof_k::Int = 0, 0
        x::Vector{Float64} = zeros(Float64, 2)
        
        for cell = 1 : ncells
            ## update FE basis evaluators
            FEBasis_∇.citem[] = cell
            update_basis!(FEBasis_∇) 
    
            ## assemble local stiffness matrix
            for j = 1 : ndofs4cell, k = j : ndofs4cell
                temp = 0
                for qp = 1 : nweights
                    temp += weights[qp] * dot(view(∇vals,:,j,qp), view(∇vals,:,k,qp))
                end
                Aloc[j,k] = temp
            end
            Aloc .*= μ * cellvolumes[cell]
    
            ## add local matrix to global matrix
            for j = 1 : ndofs4cell
                dof_j = CellDofs[j, cell]
                for k = j : ndofs4cell
                    dof_k = CellDofs[k, cell]
                    if abs(Aloc[j,k]) > 1e-15
                        # write into sparse matrix, only lines with allocations
                        rawupdateindex!(A, +, Aloc[j,k], dof_j, dof_k) 
                        if k > j
                            rawupdateindex!(A, +, Aloc[j,k], dof_k, dof_j)
                        end
                    end
                end
            end
            fill!(Aloc, 0)

            ## assemble right-hand side
            update_trafo!(L2G, cell)
            for j = 1 : ndofs4cell
                ## right-hand side
                temp = 0
                for qp = 1 : nweights
                    ## get global x for quadrature point
                    eval_trafo!(x, L2G, xref[qp])
                    ## evaluate (f(x), v_j(x))
                    temp += weights[qp] * idvals[1, j, qp] * f(x)
                end
                ## write into global vector
                dof_j = CellDofs[j, cell]
                b[dof_j] += temp * cellvolumes[cell]
            end
        end
    end
    barrier(EG, L2G)
    flush!(A)
end

assemble! (generic function with 2 methods)