Local-Global Solver API

[termi-official]

Feature List

• Local variable handler
• Custom quadrature rules via semidiscretize
• MultiLevelNewtonRaphson
• Stiff DAE solver baseline (via wrapper API)

TODOs

• Streamline multilevel setup
• Provide caches for the local solve
• Fix naming inconsistencies between local and internal
• Communication of convergence information between outer and inner solver

[termi-official]

Found some snippets. Probably part of this PR.

# Ferrite-specific
struct MaterialCache{M, DH, CV, FV}
    material::M
    dh::DH
    cv::CV
    fv::FV
end

# Standard assembly loop for global problem
function stress!(du, u, q, cache::MaterialCache, t)
    K = allocate_matrix(cache.dh)
    f = zeros(ndofs(cache.dh))
    # Allocate the element stiffness matrix
    n_basefuncs = getnbasefunctions(cache.cv)
    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(cache.dh)
        # Update the shape function gradients based on the cell coordinates
        reinit!(cache.cv, cell)
        # Reset the element stiffness matrix
        fill!(ke, 0.0)
        fill!(fe, 0.0)
        # Compute element contribution
        assemble_cell!(ke, fe, u, q, cell, cache.cv, material)
        # Assemble ke into K
        assemble!(assembler, celldofs(cell), ke, fe)
    end
    for face in FacetIterator(cache.dh, getfacetset(cache.dh.grid,"right"))
        # Update the facetvalues to the correct facet number
        reinit!(cache.fv, face)
        # Reset the temporary array for the next facet
        fill!(fe, 0.0)
        # Access the cell's coordinates
        cell_coordinates = getcoordinates(face)
        for qp in 1:getnquadpoints(cache.fv)
            # Calculate the global coordinate of the quadrature point.
            x = spatial_coordinate(cache.fv, qp, cell_coordinates)
            tₚ = cache.material.traction(x, t)
            # Get the integration weight for the current quadrature point.
            dΓ = getdetJdV(cache.fv, qp)
            for i in 1:getnbasefunctions(cache.fv)
                Nᵢ = shape_value(cache.fv, qp, i)
                fe[i] += tₚ * Nᵢ * dΓ
            end
        end
        # Add the local contributions to the correct indices in the global external force vector
        assemble!(f, celldofs(face), fe)
    end
    du .= K*u .- f
end

function stress_Ju!(K, u, q, cache::MaterialCache, t)
    # Allocate the element stiffness matrix
    n_basefuncs = getnbasefunctions(cache.cv)
    ke = zeros(n_basefuncs, n_basefuncs)
    # Create an assembler
    assembler = start_assemble(K; fillzero=false)
    # Loop over all cells
    for cell in CellIterator(dh)
        # Update the shape function gradients based on the cell coordinates
        reinit!(cache.cv, cell)
        # Reset the element stiffness matrix
        fill!(ke, 0.0)
        # Compute element contribution
        assemble_cell_Ju!(ke, u, q, cell, cache.cv, cache.material)
        # Assemble ke into K
        assemble!(assembler, celldofs(cell), ke)
    end
end

function stress_Jq!(K, u, q, cache::MaterialCache, t)
    # Allocate the element stiffness matrix
    n_basefuncs = getnbasefunctions(cache.cv)
    ke = zeros(n_basefuncs, n_basefuncs)
    # Create an assembler
    assembler = start_assemble(K)
    # Loop over all cells
    for cell in CellIterator(dh)
        # Update the shape function gradients based on the cell coordinates
        reinit!(cache.cv, cell)
        # Reset the element stiffness matrix
        fill!(ke, 0.0)
        # Compute element contribution
        assemble_cell_Jq!(ke, u, q, cell, cache.cv, cache.material)
        # Assemble ke into K
        assemble!(assembler, celldofs(cell), ke)
    end
end