text:
module Bem2d
using LaTeXStrings
using Colors
using ColorSchemes
using PyCall
using PyPlot
using Infiltrator

export Elements
mutable struct Elements
    x1::Array{Float64,1}
    y1::Array{Float64,1}
    x2::Array{Float64,1}
    y2::Array{Float64,1}
    name::Array{String,1}
    uxconst::Array{Float64,1}
    uyconst::Array{Float64,1}
    uxquad::Array{Float64,2}
    uyquad::Array{Float64,2}
    angle::Array{Float64,1}
    length::Array{Float64,1}
    halflength::Array{Float64,1}
    xcenter::Array{Float64,1}
    ycenter::Array{Float64,1}
    rotmat::Array{Float64,3}
    rotmatinv::Array{Float64,3}
    xnormal::Array{Float64,1}
    ynormal::Array{Float64,1}
    xnodes::Array{Float64,2}
    ynodes::Array{Float64,2}
    a::Array{Float64,1}
    b::Array{Float64,1}
    normalstress::Array{Float64,1}
    endidx::Int64
    Elements(maxidx) = new(fill(NaN, maxidx), # x1
                           fill(NaN, maxidx), # y1
                           fill(NaN, maxidx), # x2
                           fill(NaN, maxidx), # y2
                           fill("", maxidx), # names
                           fill(NaN, maxidx), # uxglobal::Array{Float64, 1}
                           fill(NaN, maxidx), # uyglobal::Array{Float64, 1}
                           fill(NaN, maxidx, 3), # uxglobalquad::Array{Float64, 2}
                           fill(NaN, maxidx, 3), # uyglobalquad::Array{Float64, 2}
                           fill(NaN, maxidx), # angle::Array{Float64, 1}
                           fill(NaN, maxidx), # length::Array{Float64, 1}
                           fill(NaN, maxidx), # halflength::Array{Float64, 1}
                           fill(NaN, maxidx), # xcenter::Array{Float64, 1}
                           fill(NaN, maxidx), # ycenter::Array{Float64, 1}
                           fill(NaN, maxidx, 2, 2), # rotationmatrix::Array{Float64, 3}
                           fill(NaN, maxidx, 2, 2), # rotationmatrixinverse::Array{Float64, 3}
                           fill(NaN, maxidx), # xnormal::Array{Float64, 1}
                           fill(NaN, maxidx), # ynormal::Array{Float64, 1}
                           fill(NaN, maxidx, 3), # xnodes
                           fill(NaN, maxidx, 3), # ynodes
                           fill(NaN, maxidx), # a::Array{Float64, 1}
                           fill(NaN, maxidx), # b::Array{Float64, 1}
                           fill(NaN, maxidx), # sigman::Array{Float64, 1}
                           0) # endidx
end

export discretizedarc
"""
    discretized_arc(thetastart, thetaend, radius, n_pts)

Generate regularly spaced eleemnts along an curved arc.
"""
function discretizedarc(thetastart, thetaend, radius, n_pts)
    # Create geometry of discretized arc
    thetarange = collect(LinRange(thetastart, thetaend, n_pts + 1))
    x = @. radius * cos(thetarange)
    y = @. radius * sin(thetarange)
    x1 = x[1:1:end-1]
    x2 = x[2:1:end]
    y1 = y[1:1:end-1]
    y2 = y[2:1:end]
    return x1, y1, x2, y2
end

export addelsez!
"""
    addelsez!()

Covenience function for quicly adding elements with geometry and name only
"""
function addelsez!(els, x1, y1, x2, y2, namestring)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = namestring
    end
    standardize_elements!(els)
end

export flipud
"""
    flipud()

Just an alias to the 1-linear that replicates matlab's matrix flipud
Taken from:
https://cheatsheets.quantecon.org/#manipulating-vectors-and-matrices
"""
function flipud(mat)
    mat = reverse(mat, dims = 1)
    return(mat)
end

export interleave
"""
    interleave(vec1, vec2)

Interleave values in two arrays of equal length.

# Example
```julia-repl
julia> interleave([1, 2], [3, 4])
4-element Array{Int64,1}:
 1
 3
 2
 4
```
"""
function interleave(vec1, vec2)
    return transpose([vec1 vec2])[:]
end

export meshgrid
"""
    meshgrid(xs, ys)

Replicate Matlab's meshgrid.  Makes plaid grids from two vectors.
"""
function meshgrid(xs, ys)
    xmat = repeat(xs, 1, length(ys))
    ymat = transpose(repeat(ys, 1, length(xs)))
    return xmat, ymat
end

export obsgrid
"""
    obsgrid(xmin, ymin, xmax, ymax, npts)

Create set of observation corrdinates on a regular rectangular grid.
"""
function obsgrid(xmin, ymin, xmax, ymax, npts)
    T = typeof(xmin)
    xobs, yobs = meshgrid(LinRange(T(xmin), T(xmax), npts), LinRange(T(ymin), T(ymax), npts))
    return xobs[:], yobs[:]
end

export updateendidx!
"""
    updateendidx!(els)

In place determinate of last index of structures els in which there is data
"""
function updateendidx!(els)
    els.endidx = findall(isnan, els.x1)[1] - 1
    return nothing
end

export discretizedline
"""
    discretizedline(xstart, ystart, xend, yend, nelements)

Calculated start and endpoint coordinates for a line discretized
into nelements
"""
function discretizedline(xstart, ystart, xend, yend, nelements)
    npts = nelements + 1
    x = range(xstart, stop = xend, length = npts)
    y = range(ystart, stop = yend, length = npts)
    x1 = x[1:1:end - 1]
    y1 = y[1:1:end - 1]
    x2 = x[2:1:end]
    y2 = y[2:1:end]
    return collect(x1), collect(y1), collect(x2), collect(y2)
end

export multmatvec
"""
    multmatvec(mats, vec1, vec2)

Matrix vector multiply. Multiple matrices.  Multiple vectors.
"""
function multmatvec(mats, vec1, vec2)
    newvec1 = zeros(size(vec1))
    newvec2 = zeros(size(vec2))
    for i in 1:length(vec1)
        newvec1[i], newvec2[i] = mats[i, :, :] * [vec1[i] ; vec2[i]]
    end
    return newvec1, newvec2
end

export multmatvec!
"""
    multmatvec!(mats, vec1, vec2)

In place matrix vector multiply. Multiple matrices.  Multiple vectors.
"""
function multmatvec!(vx, vy, mats, vec1, vec2)
    @inbounds for i in 1:length(vec1)
        @views vx[i] = mats[i, 1, 1] * vec1[i] + mats[i, 1, 2] * vec2[i]
        @views vy[i] = mats[i, 2, 1] * vec1[i] + mats[i, 2, 2] * vec2[i]
    end
    return nothing
end

export multmatvecquad!
"""
    multmatvecquad!(mats, vec1, vec2)

In place matrix vector multiply. Multiple matrices.  Multiple vectors.
Built for quadratic case.
"""
function multmatvecquad!(vx, vy, mats, vec1, vec2)
    for i in 1:length(vec1)
        matidx = Int64(floor((i - 1) / 3) + 1) # Change w/ every 3rd node
        @views vx[i] = mats[matidx, 1, 1] * vec1[i] + mats[matidx, 1, 2] * vec2[i]
        @views vy[i] = mats[matidx, 2, 1] * vec1[i] + mats[matidx, 2, 2] * vec2[i]
    end
    return nothing
end

export multmatsinglevec
"""
    multmatvecsingle(mats, vec1, vec2)

Matrix vector multiply. Single matrix.  Multiple vectors.
"""
function multmatsinglevec(mat, vec1, vec2)
    newvec1 = zeros(size(vec1))
    newvec2 = zeros(size(vec2))
    for i in 1:length(vec1)
        newvec1[i], newvec2[i] = mat * [vec1[i] ; vec2[i]]
    end
    return newvec1, newvec2
end

export multmatvecsingle
"""
    multmatvecsingle(mats, vec1, vec2)

Matrix vector multiply. Multiple matrices.  Single set of vectors.
"""
function multmatvecsingle(mats, vec1, vec2)
    newvec1 = zeros(size(mats)[1])
    newvec2 = zeros(size(mats)[1])
    for i in 1:size(mats)[1]
        newvec1[i], newvec2[i] = mats[i, :, :] * [vec1 ; vec2]
    end
    return newvec1, newvec2
end

export constdispstress
"""
    constdispstress(fun2dispstress, x, y, els, idx, xcomp, ycomp, mu, nu)

Calculate displacements in global reference frame for constant elements.
"""
function constdispstress(fun2dispstress, x, y, els, idx, xcomp, ycomp, mu, nu)
    disp, stress = zeros(length(x), 2), zeros(length(x), 3)
    _x, _y = zeros(length(x)), zeros(length(x))
    _xcomp, _ycomp = zeros(1), zeros(1)
    f = zeros(length(x), 7)
    _disp, _stress = zeros(length(x), 2), zeros(length(x), 3)
    @inbounds @simd for j in 1:length(idx)
        @views _x, _y = multmatsinglevec(els.rotmatinv[idx[j], :, :], x .- els.xcenter[idx[j]], y .- els.ycenter[idx[j]])
        @views _xcomp, _ycomp = els.rotmatinv[idx[j], :, :] * [xcomp[j] ; ycomp[j]]
        @views f = constkernel(_x, _y, els.halflength[idx[j]], nu)
        _disp, _stress = fun2dispstress(_xcomp, _ycomp, f, _y, mu, nu)
        @views _disp, _stress = rotdispstress(_disp, _stress, els.rotmat[idx[j], :, :])
        disp += _disp
        stress += _stress
    end
    return disp, stress
end

export quaddispstress
"""
    quaddispstress(fun2dispstress, x, y, els, idx, xcomp, ycomp, mu, nu)

Calculate displacements in global reference frame for quadratic elements.
Far-field case.
"""
function quaddispstress(fun2dispstress, x, y, els, idx, xcomp, ycomp, mu, nu)
    disp = zeros(length(x), 2)
    stress = zeros(length(x), 3)

    for j in 1:length(idx)
        _x, _y = multmatsinglevec(els.rotmatinv[idx[j], :, :], x .- els.xcenter[idx[j]], y .- els.ycenter[idx[j]])
        f = quadkernel_farfield(_x, _y, els.halflength[idx[j]], nu)
        _xcomp, _ycomp = multmatsinglevec(els.rotmatinv[idx[j], :, :], xcomp[j, :], ycomp[j, :])

        # TODO: I think I need to get xnodes in a local coordinate system
        # This should be a translate and a rotate
        # Should these matrices be moved to standardize_elements?
        xnodes, ynodes = multmatsinglevec(els.rotmatinv[idx[j], :, :], els.xnodes[idx[j], :] .- els.xcenter[idx[j]], els.ynodes[idx[j], :] .- els.ycenter[idx[j]])
        phix = slip2coef(xnodes, _xcomp, els.halflength[idx[j]])
        phiy = slip2coef(xnodes, _ycomp, els.halflength[idx[j]])
        for i in 1:3
            # _disp, _stress = fun2dispstress(_xcomp[i], _ycomp[i], f[:, :, i], _y, mu, nu)
            _disp, _stress = fun2dispstress(phix[i], phiy[i], f[:, :, i], _y, mu, nu)
            _disp, _stress = rotdispstress(_disp, _stress, els.rotmat[idx[j], :, :])
            disp += _disp
            stress += _stress
        end
    end
    return disp, stress
end

export quaddispstresscoincident
"""
    quaddispstresscoincident(fun2dispstress, x, y, els, idx, xcomp, ycomp, mu, nu)

Calculate displacements in global reference frame for quadratic elements.
Coincident case.
"""
function quaddispstresscoincident(fun2dispstress, x, y, els, idx, xcomp, ycomp, mu, nu)
    disp = zeros(length(x), 2)
    stress = zeros(length(x), 3)
    for j in 1:length(idx)
        # Rotate and translate into SC coordinate system
        _x, _y = multmatsinglevec(els.rotmatinv[idx[j], :, :], x .- els.xcenter[idx[j]], y .- els.ycenter[idx[j]])
        f = quadkernel_coincident(els.halflength[idx[j]], nu)
        _xcomp, _ycomp = multmatsinglevec(els.rotmatinv[idx[j], :, :], xcomp[j, :], ycomp[j, :])
        xnodes, ynodes = multmatsinglevec(els.rotmatinv[idx[j], :, :], els.xnodes[idx[j], :] .- els.xcenter[idx[j]], els.ynodes[idx[j], :] .- els.ycenter[idx[j]])
        phix = slip2coef(xnodes, _xcomp, els.halflength[idx[j]])
        phiy = slip2coef(xnodes, _ycomp, els.halflength[idx[j]])
        for i in 1:3
            # _disp, _stress = fun2uσ(_xcomp[i], _ycomp[i], f[:, :, i], _y, mu, nu)
            _disp, _stress = fun2dispstress(phix[i], phiy[i], f[:, :, i], _y, mu, nu)
            _disp, _stress = rotdispstress(_disp, _stress, els.rotmat[idx[j], :, :])
            disp += _disp
            stress += _stress
        end
    end
    return disp, stress
end

"""
    constkernel(x, y, a, nu)

Calculate constant slip kernels from Starfield and Crouch, pages 49 and 82.
"""
function constkernel(x, y, a, nu)
    f = zeros(length(x), 7)
    leadingconst = 1/(4*pi*(1-nu))
    for i in 1:length(x)
        # f(x) = f_1
        f[i, 1] = -leadingconst * (y[i] * (atan(y[i], (x[i]-a)) - atan(y[i], (x[i]+a))) - (x[i] - a) * log(sqrt((x[i] - a)^2 + y[i]^2)) + (x[i]+a) * log(sqrt((x[i]+a)^2 + y[i]^2)))

        # df/dy = f_2
        f[i, 2] = -leadingconst * ((atan(y[i], (x[i]-a)) - atan(y[i], (x[i]+a))))

        # df/dx = f_3
        f[i, 3] = leadingconst * (log(sqrt((x[i]-a)^2 + y[i]^2)) - log(sqrt((x[i]+a)^2 + y[i]^2)))

        # d2f/dxdy = f_4
        f[i, 4] = leadingconst * (y[i] / ((x[i]-a)^2 + y[i]^2) - y[i] / ((x[i]+a)^2 + y[i]^2))

        # d2f/dxdx = -d2f/dydy = f_5
        f[i, 5] = leadingconst * ((x[i]-a) / ((x[i]-a)^2 + y[i]^2) - (x[i]+a) / ((x[i]+a)^2 + y[i]^2))

        # d3f/dxdydy = -d3f/dxdxdx = f_6
        f[i, 6] = leadingconst * (((x[i]-a)^2 - y[i]^2) / ((x[i]-a)^2 + y[i]^2)^2 - ((x[i]+a)^2 - y[i]^2) / ((x[i]+a)^2 + y[i]^2)^2)

        # d3f/dydydy = -d3f/dxdxdy = f7
        f[i, 7] = 2 * y[i] / (4 * pi * (1 - nu)) * ((x[i]-a) / ((x[i]-a)^2 + y[i]^2)^2 - (x[i]+a) / ((x[i]+a)^2 + y[i]^2)^2)
    end
    return f
end

export quadkernel_coincident
"""
    quadkernel_coincident(a, nu)

Farfield 3-node quadratic kernels for substitution into CS kernels
f has dimensions [nobs, f=7, shapefunctions=3]

Note: These were calculated symbolically using Maxima and then
reduced to eliminate singularites by zeroing out diverging terms.
This was very tedious.
"""
function quadkernel_coincident(a, nu)
    # TODO: Should I swap order: nodeidx, f, shapefunction Probably
    f = zeros(7, 3, 3)

    f[1, 1, 1] = -5 / 144 * a * log(25 / 9 * a^2) / (pi - pi * nu) - 17 / 288 * a * log(1 / 9 * a^2) / (pi - pi * nu) + 1 / 12 * a / (pi - pi * nu)
    f[1, 2, 1] = -25 / 288 * a * log(25 / 9 * a^2) / (pi - pi * nu) + 7 / 288 * a * log(1 / 9 * a^2) / (pi - pi * nu) + 1 / 12 * a / (pi - pi * nu)
    f[1, 3, 1] = -25 / 288 * a * log(25 / 9 * a^2) / (pi - pi * nu) - 1 / 144 * a * log(1 / 9 * a^2) / (pi - pi * nu) - 1 / 6 * a / (pi - pi * nu)
    f[1, 1, 2] = -3 / 16 * a * log(a) / (pi - pi * nu) - 1 / 8 * a / (pi - pi * nu)
    f[1, 2, 2] = -1 / 8 * a * log(a) / (pi - pi * nu) + 1 / 4 * a / (pi - pi * nu)
    f[1, 3, 2] = -3 / 16 * a * log(a) / (pi - pi * nu) - 1 / 8 * a / (pi - pi * nu)
    f[1, 1, 3] = -25 / 288 * a * log(25 / 9 * a^2) / (pi - pi * nu) - 1 / 144 * a * log(1 / 9 * a^2) / (pi - pi * nu) - 1 / 6 * a / (pi - pi * nu)
    f[1, 2, 3] = -25 / 288 * a * log(25 / 9 * a^2) / (pi - pi * nu) + 7 / 288 * a * log(1 / 9 * a^2) / (pi - pi * nu) + 1 / 12 * a / (pi - pi * nu)
    f[1, 3, 3] = -5 / 144 * a * log(25 / 9 * a^2) / (pi - pi * nu) - 17 / 288 * a * log(1 / 9 * a^2) / (pi - pi * nu) + 1 / 12 * a / (pi - pi * nu)

    f[2, 1, 1] = 1 / 4 / (nu - 1)
    f[2, 2, 1] = 0
    f[2, 3, 1] = 0
    f[2, 1, 2] = 0
    f[2, 2, 2] = 1 / 4 / (nu - 1)
    f[2, 3, 2] = 0
    f[2, 1, 3] = 0
    f[2, 2, 3] = 0
    f[2, 3, 3] = 1 / 4 / (nu - 1)

    f[3, 1, 1] = 1 / 8 * log(25 / 9 * a^2) / (pi - pi * nu) - 1 / 8 * log(1 / 9 * a^2) / (pi - pi * nu) - 3 / 4 / (pi - pi * nu)
    f[3, 2, 1] = 3 / 4 / (pi - pi * nu)
    f[3, 3, 1] = 0
    f[3, 1, 2] = -3 / 8 / (pi - pi * nu)
    f[3, 2, 2] = 0
    f[3, 3, 2] = 3 / 8 / (pi - pi * nu)
    f[3, 1, 3] = 0
    f[3, 2, 3] = -3 / 4 / (pi - pi * nu)
    f[3, 3, 3] = -1 / 8 * log(25 / 9 * a^2) / (pi - pi * nu) + 1 / 8 * log(1 / 9 * a^2) / (pi - pi * nu) + 3 / 4 / (pi - pi * nu)

    f[4, 1, 1] = -9 / 16 / (a * nu - a)
    f[4, 2, 1] = 3 / 4 / (a * nu - a)
    f[4, 3, 1] = -3 / 16 / (a * nu - a)
    f[4, 1, 2] = -3 / 16 / (a * nu - a)
    f[4, 2, 2] = 0
    f[4, 3, 2] = 3 / 16 / (a * nu - a)
    f[4, 1, 3] = 3 / 16 / (a * nu - a)
    f[4, 2, 3] = -3 / 4 / (a * nu - a)
    f[4, 3, 3] = 9 / 16 / (a * nu - a)

    f[5, 1, 1] = 9 / 32 * log(25 / 9 * a^2) / (pi * a * nu - pi * a) - 9 / 32 * log(1 / 9 * a^2) / (pi * a * nu - pi * a) + 27 / 80 / (pi * a * nu - pi * a)
    f[5, 2, 1] = -3 / 8 * log(25 / 9 * a^2) / (pi * a * nu - pi * a) + 3 / 8 * log(1 / 9 * a^2) / (pi * a * nu - pi * a) + 9 / 8 / (pi * a * nu - pi * a)
    f[5, 3, 1] = 3 / 32 * log(25 / 9 * a^2) / (pi * a * nu - pi * a) - 3 / 32 * log(1 / 9 * a^2) / (pi * a * nu - pi * a) - 9 / 16 / (pi * a * nu - pi * a)
    f[5, 1, 2] = -9 / 16 / (pi * a * nu - pi * a)
    f[5, 2, 2] = 13 / 8 / (pi * a * nu - pi * a)
    f[5, 3, 2] = -9 / 16 / (pi * a * nu - pi * a)
    f[5, 1, 3] = 3 / 32 * log(25 / 9 * a^2) / (pi * a * nu - pi * a) - 3 / 32 * log(1 / 9 * a^2) / (pi * a * nu - pi * a) - 9 / 16 / (pi * a * nu - pi * a)
    f[5, 2, 3] = -3 / 8 * log(25 / 9 * a^2) / (pi * a * nu - pi * a) + 3 / 8 * log(1 / 9 * a^2) / (pi * a * nu - pi * a) + 9 / 8 / (pi * a * nu - pi * a)
    f[5, 3, 3] = 9 / 32 * log(25 / 9 * a^2) / (pi * a * nu - pi * a) - 9 / 32 * log(1 / 9 * a^2) / (pi * a * nu - pi * a) + 27 / 80 / (pi * a * nu - pi * a)

    f[6, 1, 1] = 9 / 32 * log(25 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) - 9 / 32 * log(1 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) + 621 / 100 / (pi * a^2 * nu - pi * a^2)
    f[6, 2, 1] = -9 / 16 * log(25 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) + 9 / 16 * log(1 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) - 27 / 5 / (pi * a^2 * nu - pi * a^2)
    f[6, 3, 1] = 9 / 32 * log(25 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) - 9 / 32 * log(1 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) + 27 / 20 / (pi * a^2 * nu - pi * a^2)
    f[6, 1, 2] = 3 / 4 / (pi * a^2 * nu - pi * a^2)
    f[6, 2, 2] = 0
    f[6, 3, 2] = -3 / 4 / (pi * a^2 * nu - pi * a^2)
    f[6, 1, 3] = -9 / 32 * log(25 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) + 9 / 32 * log(1 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) - 27 / 20 / (pi * a^2 * nu - pi * a^2)
    f[6, 2, 3] = 9 / 16 * log(25 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) - 9 / 16 * log(1 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) + 27 / 5 / (pi * a^2 * nu - pi * a^2)
    f[6, 3, 3] = -9 / 32 * log(25 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) + 9 / 32 * log(1 / 9 * a^2) / (pi * a^2 * nu - pi * a^2) - 621 / 100 / (pi * a^2 * nu - pi * a^2)

    f[7, 1, 1] = -9 / 16 / (a^2 * nu - a^2)
    f[7, 2, 1] = 9 / 8 / (a^2 * nu - a^2)
    f[7, 3, 1] = -9 / 16 / (a^2 * nu - a^2)
    f[7, 1, 2] = -9 / 16 / (a^2 * nu - a^2)
    f[7, 2, 2] = 9 / 8 / (a^2 * nu - a^2)
    f[7, 3, 2] = -9 / 16 / (a^2 * nu - a^2)
    f[7, 1, 3] = -9 / 16 / (a^2 * nu - a^2)
    f[7, 2, 3] = 9 / 8 / (a^2 * nu - a^2)
    f[7, 3, 3] = -9 / 16 / (a^2 * nu - a^2)

    f = permutedims(f, [3, 1, 2]) # This is just a hack until I get around to putting in the correct indices above
    return f
end

export quadkernel_farfield
"""
    quadkernel_farfield(x, y, a, nu)

Farfield 3-node quadratic kernels for substitution into CS kernels
f has dimensions [nobs, f=7, shapefunctions=3]

Note: These were calculated symbolically using Maxima.
"""
function quadkernel_farfield(x, y, a, nu)
    f = zeros(length(x), 7, 3)
    for i in 1:length(x)
        # term1 = atan((a - x[i]) / y[i])
        # term2 = atan((a + x[i]) / y[i])
        term1 = pi / 2 * sign(y[i] / (a - x[i])) - atan(y[i] / (a - x[i]))
        term2 = pi / 2 * sign(y[i] / (a + x[i])) - atan(y[i] / (a + x[i]))

        f[i, 1, 1] = -1 / 64 * ( 6 * y[i]^3 * (term2 + term1) - 6 * a^3 * log(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2) - 8 * a^3 - 12 * a^2 * x[i] + 12 * a * x[i]^2 - 12 * a * y[i]^2 + 6 * ( (2 * a * x[i] - 3 * x[i]^2) * term2 + (2 * a * x[i] - 3 * x[i]^2) * term1 ) * y[i] + 3 * (a * x[i]^2 - x[i]^3 - (a - 3 * x[i]) * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) - 3 * (a * x[i]^2 - x[i]^3 - (a - 3 * x[i]) * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2)) ) / (pi * a^2 * nu - pi * a^2)
        f[i, 1, 2] = 1 / 32 * (6 * y[i]^3 * (term2 + term1) + a^3 * log(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2) + a^3 * log(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2) - 8 * a^3 + 12 * a * x[i]^2 - 12 * a * y[i]^2 + 2 * ((4 * a^2 - 9 * x[i]^2) * term2 + (4 * a^2 - 9 * x[i]^2) * term1) * y[i] + (4 * a^2 * x[i] - 3 * x[i]^3 + 9 * x[i] * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) - (4 * a^2 * x[i] - 3 * x[i]^3 + 9 * x[i] * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^2 * nu - pi * a^2)
        f[i, 1, 3] = -1 / 64 * (6 * y[i]^3 * (term2 + term1) - 6 * a^3 * log(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2) - 8 * a^3 + 12 * a^2 * x[i] + 12 * a * x[i]^2 - 12 * a * y[i]^2 - 6 * ((2 * a * x[i] + 3 * x[i]^2) * term2 + (2 * a * x[i] + 3 * x[i]^2) * term1) * y[i] - 3 * (a * x[i]^2 + x[i]^3 - (a + 3 * x[i]) * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) + 3 * (a * x[i]^2 + x[i]^3 - (a + 3 * x[i]) * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^2 * nu - pi * a^2)

        f[i, 2, 1] = -3 / 32 * (3 * y[i]^2 * (term2 + term1) - (a - 3 * x[i]) * y[i] * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) + (a - 3 * x[i]) * y[i] * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2)) - 6 * a * y[i] + (2 * a * x[i] - 3 * x[i]^2) * term2 + (2 * a * x[i] - 3 * x[i]^2) * term1) / (pi * a^2 * nu - pi * a^2)
        f[i, 2, 2] = 1 / 16 * (9 * y[i]^2 * (term2 + term1) + 9 * x[i] * y[i] * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) - 9 * x[i] * y[i] * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2)) - 18 * a * y[i] + (4 * a^2 - 9 * x[i]^2) * term2 + (4 * a^2 - 9 * x[i]^2) * term1) / (pi * a^2 * nu - pi * a^2)
        f[i, 2, 3] = -3 / 32 * (3 * y[i]^2 * (term2 + term1) + (a + 3 * x[i]) * y[i] * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) - (a + 3 * x[i]) * y[i] * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2)) - 6 * a * y[i] - (2 * a * x[i] + 3 * x[i]^2) * term2 - (2 * a * x[i] + 3 * x[i]^2) * term1) / (pi * a^2 * nu - pi * a^2)

        f[i, 3, 1] = 3 / 64 * (8 * a^2 - 12 * a * x[i] - 4 * ((a - 3 * x[i]) * term2 + (a - 3 * x[i]) * term1) * y[i] - (2 * a * x[i] - 3 * x[i]^2 + 3 * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) + (2 * a * x[i] - 3 * x[i]^2 + 3 * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^2 * nu - pi * a^2)
        f[i, 3, 2] = 1 / 32 * (36 * a * x[i] - 36 * (x[i] * term2 + x[i] * term1) * y[i] + (4 * a^2 - 9 * x[i]^2 + 9 * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) - (4 * a^2 - 9 * x[i]^2 + 9 * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^2 * nu - pi * a^2)
        f[i, 3, 3] = -3 / 64 * (8 * a^2 + 12 * a * x[i] - 4 * ((a + 3 * x[i]) * term2 + (a + 3 * x[i]) * term1) * y[i] - (2 * a * x[i] + 3 * x[i]^2 - 3 * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) + (2 * a * x[i] + 3 * x[i]^2 - 3 * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^2 * nu - pi * a^2)

        f[i, 4, 1] = 3 / 32 * (4 * a^2 * y[i]^3 - 2 * ((a - 3 * x[i]) * term2 + (a - 3 * x[i]) * term1) * y[i]^4 - 4 * ((a^3 - 3 * a^2 * x[i] + a * x[i]^2 - 3 * x[i]^3) * term2 + (a^3 - 3 * a^2 * x[i] + a * x[i]^2 - 3 * x[i]^3) * term1) * y[i]^2 + 4 * (a^4 - 3 * a^3 * x[i] + a^2 * x[i]^2) * y[i] - 2 * (a^5 - 3 * a^4 * x[i] - 2 * a^3 * x[i]^2 + 6 * a^2 * x[i]^3 + a * x[i]^4 - 3 * x[i]^5) * term2 - 2 * (a^5 - 3 * a^4 * x[i] - 2 * a^3 * x[i]^2 + 6 * a^2 * x[i]^3 + a * x[i]^4 - 3 * x[i]^5) * term1 - 3 * (y[i]^5 + 2 * (a^2 + x[i]^2) * y[i]^3 + (a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * y[i]) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) + 3 * (y[i]^5 + 2 * (a^2 + x[i]^2) * y[i]^3 + (a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * y[i]) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^6 * nu - pi * a^6 + (pi * a^2 * nu - pi * a^2) * x[i]^4 + (pi * a^2 * nu - pi * a^2) * y[i]^4 - 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2 + 2 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^2)
        f[i, 4, 2] = 1 / 16 * (20 * a^3 * x[i] * y[i] - 18 * (x[i] * term2 + x[i] * term1) * y[i]^4 - 36 * ((a^2 * x[i] + x[i]^3) * term2 + (a^2 * x[i] + x[i]^3) * term1) * y[i]^2 - 18 * (a^4 * x[i] - 2 * a^2 * x[i]^3 + x[i]^5) * term2 - 18 * (a^4 * x[i] - 2 * a^2 * x[i]^3 + x[i]^5) * term1 + 9 * (y[i]^5 + 2 * (a^2 + x[i]^2) * y[i]^3 + (a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * y[i]) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) - 9 * (y[i]^5 + 2 * (a^2 + x[i]^2) * y[i]^3 + (a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * y[i]) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^6 * nu - pi * a^6 + (pi * a^2 * nu - pi * a^2) * x[i]^4 + (pi * a^2 * nu - pi * a^2) * y[i]^4 - 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2 + 2 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^2)
        f[i, 4, 3] = -3 / 32 * (4 * a^2 * y[i]^3 - 2 * ((a + 3 * x[i]) * term2 + (a + 3 * x[i]) * term1) * y[i]^4 - 4 * ((a^3 + 3 * a^2 * x[i] + a * x[i]^2 + 3 * x[i]^3) * term2 + (a^3 + 3 * a^2 * x[i] + a * x[i]^2 + 3 * x[i]^3) * term1) * y[i]^2 + 4 * (a^4 + 3 * a^3 * x[i] + a^2 * x[i]^2) * y[i] - 2 * (a^5 + 3 * a^4 * x[i] - 2 * a^3 * x[i]^2 - 6 * a^2 * x[i]^3 + a * x[i]^4 + 3 * x[i]^5) * term2 - 2 * (a^5 + 3 * a^4 * x[i] - 2 * a^3 * x[i]^2 - 6 * a^2 * x[i]^3 + a * x[i]^4 + 3 * x[i]^5) * term1 + 3 * (y[i]^5 + 2 * (a^2 + x[i]^2) * y[i]^3 + (a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * y[i]) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) - 3 * (y[i]^5 + 2 * (a^2 + x[i]^2) * y[i]^3 + (a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * y[i]) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^6 * nu - pi * a^6 + (pi * a^2 * nu - pi * a^2) * x[i]^4 + (pi * a^2 * nu - pi * a^2) * y[i]^4 - 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2 + 2 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^2)

        f[i, 5, 1] = 3 / 32 * (6 * y[i]^5 * (term2 + term1) - 6 * a^5 - 4 * a^4 * x[i] + 18 * a^3 * x[i]^2 + 4 * a^2 * x[i]^3 - 12 * a * x[i]^4 - 12 * a * y[i]^4 + 12 * ((a^2 + x[i]^2) * term2 + (a^2 + x[i]^2) * term1) * y[i]^3 - 2 * (9 * a^3 - 2 * a^2 * x[i] + 12 * a * x[i]^2) * y[i]^2 + 6 * ((a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * term2 + (a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * term1) * y[i] - (a^5 - 3 * a^4 * x[i] - 2 * a^3 * x[i]^2 + 6 * a^2 * x[i]^3 + a * x[i]^4 - 3 * x[i]^5 + (a - 3 * x[i]) * y[i]^4 + 2 * (a^3 - 3 * a^2 * x[i] + a * x[i]^2 - 3 * x[i]^3) * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) + (a^5 - 3 * a^4 * x[i] - 2 * a^3 * x[i]^2 + 6 * a^2 * x[i]^3 + a * x[i]^4 - 3 * x[i]^5 + (a - 3 * x[i]) * y[i]^4 + 2 * (a^3 - 3 * a^2 * x[i] + a * x[i]^2 - 3 * x[i]^3) * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^6 * nu - pi * a^6 + (pi * a^2 * nu - pi * a^2) * x[i]^4 + (pi * a^2 * nu - pi * a^2) * y[i]^4 - 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2 + 2 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^2)
        f[i, 5, 2] = -1 / 16 * (18 * y[i]^5 * (term2 + term1) - 26 * a^5 + 62 * a^3 * x[i]^2 - 36 * a * x[i]^4 - 36 * a * y[i]^4 + 36 * ((a^2 + x[i]^2) * term2 + (a^2 + x[i]^2) * term1) * y[i]^3 - 2 * (31 * a^3 + 36 * a * x[i]^2) * y[i]^2 + 18 * ((a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * term2 + (a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * term1) * y[i] + 9 * (a^4 * x[i] - 2 * a^2 * x[i]^3 + x[i]^5 + x[i] * y[i]^4 + 2 * (a^2 * x[i] + x[i]^3) * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) - 9 * (a^4 * x[i] - 2 * a^2 * x[i]^3 + x[i]^5 + x[i] * y[i]^4 + 2 * (a^2 * x[i] + x[i]^3) * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^6 * nu - pi * a^6 + (pi * a^2 * nu - pi * a^2) * x[i]^4 + (pi * a^2 * nu - pi * a^2) * y[i]^4 - 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2 + 2 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^2)
        f[i, 5, 3] = 3 / 32 * (6 * y[i]^5 * (term2 + term1) - 6 * a^5 + 4 * a^4 * x[i] + 18 * a^3 * x[i]^2 - 4 * a^2 * x[i]^3 - 12 * a * x[i]^4 - 12 * a * y[i]^4 + 12 * ((a^2 + x[i]^2) * term2 + (a^2 + x[i]^2) * term1) * y[i]^3 - 2 * (9 * a^3 + 2 * a^2 * x[i] + 12 * a * x[i]^2) * y[i]^2 + 6 * ((a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * term2 + (a^4 - 2 * a^2 * x[i]^2 + x[i]^4) * term1) * y[i] + (a^5 + 3 * a^4 * x[i] - 2 * a^3 * x[i]^2 - 6 * a^2 * x[i]^3 + a * x[i]^4 + 3 * x[i]^5 + (a + 3 * x[i]) * y[i]^4 + 2 * (a^3 + 3 * a^2 * x[i] + a * x[i]^2 + 3 * x[i]^3) * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) - (a^5 + 3 * a^4 * x[i] - 2 * a^3 * x[i]^2 - 6 * a^2 * x[i]^3 + a * x[i]^4 + 3 * x[i]^5 + (a + 3 * x[i]) * y[i]^4 + 2 * (a^3 + 3 * a^2 * x[i] + a * x[i]^2 + 3 * x[i]^3) * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^6 * nu - pi * a^6 + (pi * a^2 * nu - pi * a^2) * x[i]^4 + (pi * a^2 * nu - pi * a^2) * y[i]^4 - 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2 + 2 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^2)

        f[i, 6, 1] = 3 / 32 * (8 * a^8 - 24 * a^7 * x[i] - 16 * a^6 * x[i]^2 + 60 * a^5 * x[i]^3 + 8 * a^4 * x[i]^4 - 48 * a^3 * x[i]^5 + 12 * a * x[i]^7 + 12 * a * x[i] * y[i]^6 + 4 * (2 * a^4 + 12 * a^3 * x[i] + 9 * a * x[i]^3) * y[i]^4 + 4 * (4 * a^6 + 3 * a^5 * x[i] - 12 * a^4 * x[i]^2 + 9 * a * x[i]^5) * y[i]^2 - 3 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8 + y[i]^8 + 4 * (a^2 + x[i]^2) * y[i]^6 + 2 * (3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * y[i]^4 + 4 * (a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) + 3 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8 + y[i]^8 + 4 * (a^2 + x[i]^2) * y[i]^6 + 2 * (3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * y[i]^4 + 4 * (a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^10 * nu - pi * a^10 + (pi * a^2 * nu - pi * a^2) * x[i]^8 + (pi * a^2 * nu - pi * a^2) * y[i]^8 - 4 * (pi * a^4 * nu - pi * a^4) * x[i]^6 + 4 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^6 + 6 * (pi * a^6 * nu - pi * a^6) * x[i]^4 + 2 * (3 * pi * a^6 * nu - 3 * pi * a^6 + 3 * (pi * a^2 * nu - pi * a^2) * x[i]^4 + 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2) * y[i]^4 - 4 * (pi * a^8 * nu - pi * a^8) * x[i]^2 + 4 * (pi * a^8 * nu - pi * a^8 + (pi * a^2 * nu - pi * a^2) * x[i]^6 - (pi * a^4 * nu - pi * a^4) * x[i]^4 - (pi * a^6 * nu - pi * a^6) * x[i]^2) * y[i]^2)
        f[i, 6, 2] = 1 / 16 * (56 * a^7 * x[i] - 148 * a^5 * x[i]^3 + 128 * a^3 * x[i]^5 - 36 * a * x[i]^7 - 36 * a * x[i] * y[i]^6 - 12 * (8 * a^3 * x[i] + 9 * a * x[i]^3) * y[i]^4 - 4 * (a^5 * x[i] - 8 * a^3 * x[i]^3 + 27 * a * x[i]^5) * y[i]^2 + 9 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8 + y[i]^8 + 4 * (a^2 + x[i]^2) * y[i]^6 + 2 * (3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * y[i]^4 + 4 * (a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) - 9 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8 + y[i]^8 + 4 * (a^2 + x[i]^2) * y[i]^6 + 2 * (3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * y[i]^4 + 4 * (a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^10 * nu - pi * a^10 + (pi * a^2 * nu - pi * a^2) * x[i]^8 + (pi * a^2 * nu - pi * a^2) * y[i]^8 - 4 * (pi * a^4 * nu - pi * a^4) * x[i]^6 + 4 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^6 + 6 * (pi * a^6 * nu - pi * a^6) * x[i]^4 + 2 * (3 * pi * a^6 * nu - 3 * pi * a^6 + 3 * (pi * a^2 * nu - pi * a^2) * x[i]^4 + 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2) * y[i]^4 - 4 * (pi * a^8 * nu - pi * a^8) * x[i]^2 + 4 * (pi * a^8 * nu - pi * a^8 + (pi * a^2 * nu - pi * a^2) * x[i]^6 - (pi * a^4 * nu - pi * a^4) * x[i]^4 - (pi * a^6 * nu - pi * a^6) * x[i]^2) * y[i]^2)
        f[i, 6, 3] = -3 / 32 * (8 * a^8 + 24 * a^7 * x[i] - 16 * a^6 * x[i]^2 - 60 * a^5 * x[i]^3 + 8 * a^4 * x[i]^4 + 48 * a^3 * x[i]^5 - 12 * a * x[i]^7 - 12 * a * x[i] * y[i]^6 + 4 * (2 * a^4 - 12 * a^3 * x[i] - 9 * a * x[i]^3) * y[i]^4 + 4 * (4 * a^6 - 3 * a^5 * x[i] - 12 * a^4 * x[i]^2 - 9 * a * x[i]^5) * y[i]^2 + 3 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8 + y[i]^8 + 4 * (a^2 + x[i]^2) * y[i]^6 + 2 * (3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * y[i]^4 + 4 * (a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * y[i]^2) * log(abs(a^2 + 2 * a * x[i] + x[i]^2 + y[i]^2)) - 3 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8 + y[i]^8 + 4 * (a^2 + x[i]^2) * y[i]^6 + 2 * (3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * y[i]^4 + 4 * (a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * y[i]^2) * log(abs(a^2 - 2 * a * x[i] + x[i]^2 + y[i]^2))) / (pi * a^10 * nu - pi * a^10 + (pi * a^2 * nu - pi * a^2) * x[i]^8 + (pi * a^2 * nu - pi * a^2) * y[i]^8 - 4 * (pi * a^4 * nu - pi * a^4) * x[i]^6 + 4 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^6 + 6 * (pi * a^6 * nu - pi * a^6) * x[i]^4 + 2 * (3 * pi * a^6 * nu - 3 * pi * a^6 + 3 * (pi * a^2 * nu - pi * a^2) * x[i]^4 + 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2) * y[i]^4 - 4 * (pi * a^8 * nu - pi * a^8) * x[i]^2 + 4 * (pi * a^8 * nu - pi * a^8 + (pi * a^2 * nu - pi * a^2) * x[i]^6 - (pi * a^4 * nu - pi * a^4) * x[i]^4 - (pi * a^6 * nu - pi * a^6) * x[i]^2) * y[i]^2)

        f[i, 7, 1] = -3 / 16 * (3 * y[i]^8 * (term2 + term1) - 6 * a * y[i]^7 + 12 * ((a^2 + x[i]^2) * term2 + (a^2 + x[i]^2) * term1) * y[i]^6 - 6 * (4 * a^3 + 3 * a * x[i]^2) * y[i]^5 + 6 * ((3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * term2 + (3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * term1) * y[i]^4 - 2 * (15 * a^5 - 8 * a^4 * x[i] + 9 * a * x[i]^4) * y[i]^3 + 12 * ((a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * term2 + (a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * term1) * y[i]^2 - 2 * (6 * a^7 - 8 * a^6 * x[i] + 3 * a^5 * x[i]^2 + 8 * a^4 * x[i]^3 - 12 * a^3 * x[i]^4 + 3 * a * x[i]^6) * y[i] + 3 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8) * term2 + 3 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8) * term1) / (pi * a^10 * nu - pi * a^10 + (pi * a^2 * nu - pi * a^2) * x[i]^8 + (pi * a^2 * nu - pi * a^2) * y[i]^8 - 4 * (pi * a^4 * nu - pi * a^4) * x[i]^6 + 4 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^6 + 6 * (pi * a^6 * nu - pi * a^6) * x[i]^4 + 2 * (3 * pi * a^6 * nu - 3 * pi * a^6 + 3 * (pi * a^2 * nu - pi * a^2) * x[i]^4 + 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2) * y[i]^4 - 4 * (pi * a^8 * nu - pi * a^8) * x[i]^2 + 4 * (pi * a^8 * nu - pi * a^8 + (pi * a^2 * nu - pi * a^2) * x[i]^6 - (pi * a^4 * nu - pi * a^4) * x[i]^4 - (pi * a^6 * nu - pi * a^6) * x[i]^2) * y[i]^2)
        f[i, 7, 2] = 1 / 8 * (9 * y[i]^8 * (term2 + term1) - 18 * a * y[i]^7 + 36 * ((a^2 + x[i]^2) * term2 + (a^2 + x[i]^2) * term1) * y[i]^6 - 2 * (32 * a^3 + 27 * a * x[i]^2) * y[i]^5 + 18 * ((3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * term2 + (3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * term1) * y[i]^4 - 2 * (37 * a^5 + 8 * a^3 * x[i]^2 + 27 * a * x[i]^4) * y[i]^3 + 36 * ((a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * term2 + (a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * term1) * y[i]^2 - 2 * (14 * a^7 + a^5 * x[i]^2 - 24 * a^3 * x[i]^4 + 9 * a * x[i]^6) * y[i] + 9 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8) * term2 + 9 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8) * term1) / (pi * a^10 * nu - pi * a^10 + (pi * a^2 * nu - pi * a^2) * x[i]^8 + (pi * a^2 * nu - pi * a^2) * y[i]^8 - 4 * (pi * a^4 * nu - pi * a^4) * x[i]^6 + 4 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^6 + 6 * (pi * a^6 * nu - pi * a^6) * x[i]^4 + 2 * (3 * pi * a^6 * nu - 3 * pi * a^6 + 3 * (pi * a^2 * nu - pi * a^2) * x[i]^4 + 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2) * y[i]^4 - 4 * (pi * a^8 * nu - pi * a^8) * x[i]^2 + 4 * (pi * a^8 * nu - pi * a^8 + (pi * a^2 * nu - pi * a^2) * x[i]^6 - (pi * a^4 * nu - pi * a^4) * x[i]^4 - (pi * a^6 * nu - pi * a^6) * x[i]^2) * y[i]^2)
        f[i, 7, 3] = -3 / 16 * (3 * y[i]^8 * (term2 + term1) - 6 * a * y[i]^7 + 12 * ((a^2 + x[i]^2) * term2 + (a^2 + x[i]^2) * term1) * y[i]^6 - 6 * (4 * a^3 + 3 * a * x[i]^2) * y[i]^5 + 6 * ((3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * term2 + (3 * a^4 + 2 * a^2 * x[i]^2 + 3 * x[i]^4) * term1) * y[i]^4 - 2 * (15 * a^5 + 8 * a^4 * x[i] + 9 * a * x[i]^4) * y[i]^3 + 12 * ((a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * term2 + (a^6 - a^4 * x[i]^2 - a^2 * x[i]^4 + x[i]^6) * term1) * y[i]^2 - 2 * (6 * a^7 + 8 * a^6 * x[i] + 3 * a^5 * x[i]^2 - 8 * a^4 * x[i]^3 - 12 * a^3 * x[i]^4 + 3 * a * x[i]^6) * y[i] + 3 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8) * term2 + 3 * (a^8 - 4 * a^6 * x[i]^2 + 6 * a^4 * x[i]^4 - 4 * a^2 * x[i]^6 + x[i]^8) * term1) / (pi * a^10 * nu - pi * a^10 + (pi * a^2 * nu - pi * a^2) * x[i]^8 + (pi * a^2 * nu - pi * a^2) * y[i]^8 - 4 * (pi * a^4 * nu - pi * a^4) * x[i]^6 + 4 * (pi * a^4 * nu - pi * a^4 + (pi * a^2 * nu - pi * a^2) * x[i]^2) * y[i]^6 + 6 * (pi * a^6 * nu - pi * a^6) * x[i]^4 + 2 * (3 * pi * a^6 * nu - 3 * pi * a^6 + 3 * (pi * a^2 * nu - pi * a^2) * x[i]^4 + 2 * (pi * a^4 * nu - pi * a^4) * x[i]^2) * y[i]^4 - 4 * (pi * a^8 * nu - pi * a^8) * x[i]^2 + 4 * (pi * a^8 * nu - pi * a^8 + (pi * a^2 * nu - pi * a^2) * x[i]^6 - (pi * a^4 * nu - pi * a^4) * x[i]^4 - (pi * a^6 * nu - pi * a^6) * x[i]^2) * y[i]^2)
    end
    return f
end

export trac2dispstress
"""
    trac2dispstress(xcomp, ycomp, f, y, mu, nu)

Calculate interior displacements and stresses from boundary tractions
"""
function trac2dispstress(xcomp, ycomp, f, y, mu, nu)
    disp = zeros(length(y), 2)
    stress = zeros(length(y), 3)
    # _xcomp, _ycomp = -xcomp, -ycomp # For Okada consistency...causes problems...don't do this.
    _xcomp, _ycomp = xcomp, ycomp # No sign flip.  This provides a different standard than Okada but makes more sense for BEM
    # TODO: Check sigma_yy and sigma_xy for forcing in the x-direction
    for i in 1:length(y)
        disp[i, 1] = _xcomp / (2.0 * mu) * ((3.0 - 4.0 * nu) * f[i, 1] + y[i] * f[i, 2]) + _ycomp / (2.0 * mu) * (-y[i] * f[i, 3])
        disp[i, 2] = _xcomp / (2.0 * mu) * (-y[i] * f[i, 3]) + _ycomp / (2.0 * mu) * ((3.0 - 4.0 * nu) * f[i, 1] - y[i] * f[i, 2])
        stress[i, 1] = _xcomp * ((3.0 - 2.0 * nu) * f[i, 3] + y[i] * f[i, 4]) + _ycomp * (2.0 * nu * f[i, 2] + y[i] * -f[i, 5])
        # stress[i, 2] = _xcomp * (-1.0 * (1.0 - 2.0 * nu) * f[i, 3] + y[i] * f[i, 4]) + _ycomp * (2.0 * (1.0 - nu) * f[i, 2] - y[i] * -f[i, 5])
        # stress[i, 3] = _xcomp * (2.0 * (1.0 - nu) * f[i, 2] + y[i] * f[i, 5]) + _ycomp * ((1.0 - 2.0 * nu) * f[i, 3] - y[i] * f[i, 4])
        stress[i, 2] = _xcomp * (-1.0 * (1.0 - 2.0 * nu) * f[i, 3] - y[i] * f[i, 4]) + _ycomp * (2.0 * (1.0 - nu) * f[i, 2] - y[i] * -f[i, 5])
        stress[i, 3] = _xcomp * (2.0 * (1.0 - nu) * f[i, 2] + y[i] * -f[i, 5]) + _ycomp * ((1.0 - 2.0 * nu) * f[i, 3] - y[i] * f[i, 4])
    end
    return disp, stress
end

export slip2dispstress
"""
    slip2dispstress(xcomp, ycomp, f, y, mu, nu)

Calculate interior displacements and stresses from boundary slip/displacement
"""
function slip2dispstress(xcomp, ycomp, f, y, mu, nu)
    disp, stress = zeros(length(y), 2), zeros(length(y), 3)
    _xcomp, _ycomp = -xcomp, -ycomp # For Okada consistency
    # _xcomp, _ycomp = xcomp, ycomp # TODO: Neglecting the sign flip leads to errors with ...super strange
    for i in 1:length(y)
        # CS equations (5.2.3), page 81
        disp[i, 1] = _xcomp * (2.0 * (1.0 - nu) * f[i, 2] - y[i] * f[i, 5]) +
                     _ycomp * (-(1.0 - 2.0 * nu) * f[i, 3] - y[i] * f[i, 4])
        disp[i, 2] = _xcomp * ((1.0 - 2.0 * nu) * f[i, 3] - y[i] * f[i, 4]) +
                     _ycomp * (2.0 * (1 - nu) * f[i, 2] - y[i] * -f[i, 5])

        # CS equations (5.2.4), page 81
        stress[i, 1] = 2.0 * _xcomp * mu * (2.0 * f[i, 4] + y[i] * f[i, 6]) +
                       2.0 * _ycomp * mu * (-f[i, 5] + y[i] * f[i, 7])
        stress[i, 2] = 2.0 * _xcomp * mu * (-y[i] * f[i, 6]) +
                       2.0 * _ycomp * mu * (-f[i, 5] - y[i] * f[i, 7])
        stress[i, 3] = 2.0 * _xcomp * mu * (-f[i, 5] + y[i] * f[i, 7]) +
                       2.0 * _ycomp * mu * (-y[i] * f[i, 6])
    end
    return disp, stress
end

export stress2trac
"""
    stress2trac(stress, nvec)

Calculate traction from stresses and normal vector.
"""
function stress2trac(stress, nvec)
    return [stress[1] stress[3] ; stress[3] stress[2]] * nvec
end

export standardize_elements!
"""
    standardize_elements!(els)

In place updating els structure to determin derived parameters
(e.g., midpoint, length, etc.)

Note: Should be called after declaring elements and before any
BEM calculations (boundary or interior)
"""
function standardize_elements!(els)
    updateendidx!(els)
    for i in 1:els.endidx
        dx = els.x2[i] - els.x1[i]
        dy = els.y2[i] - els.y1[i]
        magnitude = sqrt(dx^2 + dy^2)
        els.angle[i] = atan(dy, dx)
        els.length[i] = magnitude
        els.halflength[i] = 0.5 * els.length[i]
        els.xcenter[i] = 0.5 * (els.x2[i] + els.x1[i])
        els.ycenter[i] = 0.5 * (els.y2[i] + els.y1[i])
        els.rotmat[i, :, :] = [cos(els.angle[i]) -sin(els.angle[i]) ; sin(els.angle[i]) cos(els.angle[i])]
        els.rotmatinv[i, :, :] = [cos(-els.angle[i]) -sin(-els.angle[i]) ; sin(-els.angle[i]) cos(-els.angle[i])]
        els.xnormal[i] = dy / magnitude
        els.ynormal[i] = -dx / magnitude
        els.xnodes[i, :] = [els.xcenter[i] - (2 / 3 * dx / 2), els.xcenter[i], els.xcenter[i] + (2 / 3 * dx / 2)]
        els.ynodes[i, :] = [els.ycenter[i] - (2 / 3 * dy / 2), els.ycenter[i], els.ycenter[i] + (2 / 3 * dy / 2)]
    end
    return nothing
end

export rotdispstress
"""
    rotdispstress(disp, stress, rotmat)

Rotate displacement vectors and stress tensors.

Note: Consider moving to StaticArrays
"""
function rotdispstress(disp, stress, rotmat)
    # TODO: If this is slow hand expand the matrix vector mulltiplies
    # TODO: inplace for speed.  Some benchmarks suggest 50x speedup!
    _disp = zeros(size(disp))
    _stress = zeros(size(stress))
    for i in 1:size(stress)[1]
        _disp[i, 1], _disp[i, 2] = rotmat * [disp[i, 1] ; disp[i, 2]]
        stresstensor = [stress[i, 1] stress[i, 3] ; stress[i, 3] stress[i, 2]]
        _stress[i, 1], _stress[i, 3], _, _stress[i, 2] = rotmat * stresstensor * transpose(rotmat)
    end
    return _disp, _stress
end

export plotelements
"""
    plotelements(els)

Plot all non-NaN parts of elements as straing lines.
"""
function plotelements(els)
    for i in 1:els.endidx
        plot([els.x1[i], els.x2[i]], [els.y1[i], els.y2[i]], "-k", linewidth = 2.0)
    end
    return nothing
end

"""
    stylesubplots(xlim, ylim)

Common styling for subplots generated by plotfields.
"""
function stylesubplots(xlim, ylim)
    gca().set_aspect("equal")
    gca().set_xlim([xlim[1], xlim[2]])
    gca().set_ylim([ylim[1], ylim[2]])
    gca().set_xticks([xlim[1], xlim[2]])
    gca().set_yticks([ylim[1], ylim[2]])
    return nothing
end

"""
    plotfields_contours(els, xobs, yobs, idx, field, title)

Contouring of a single field for six-panel plotfields.
"""
function plotfields_contours(els, xobs, yobs, idx, field, title)
    ncontours = 40
    xlim = [minimum(xobs) maximum(xobs)]
    ylim = [minimum(yobs) maximum(yobs)]
    subplot(2, 3, idx)
    scale = 1.0
    fieldmax = maximum(@.abs(field))
    contourf(xobs, yobs, reshape(field, size(xobs)), ncontours,
        vmin = -scale * fieldmax, vmax = scale * fieldmax, cmap = rycroftcmap())
    clim(-scale * fieldmax, scale * fieldmax)
    colorbar(fraction = 0.020, pad = 0.05, extend = "both")
    contour(xobs, yobs, reshape(field, size(xobs)), ncontours,
        vmin = -scale * fieldmax, vmax = scale * fieldmax, linewidths = 0.25, colors = "k")
    PyPlot.title(title)
    stylesubplots(xlim, ylim)
    plotelements(els)
    return nothing
end

export plotfields
"""
    plotfields(els, xobs, yobs, disp, stress, suptitle)

Set up six panel plots for displacement and stress fields.
"""
function plotfields(els, xobs, yobs, disp, stress, suptitle)
    figure(figsize = (25, 15))
    subplot(2, 3, 1)
    quiver(xobs[:], yobs[:], disp[:, 1], disp[:, 2], units = "width", color = "b")
    stylesubplots([minimum(xobs) maximum(xobs)], [minimum(yobs) maximum(yobs)])
    plotelements(els)
    PyPlot.title("displacements")
    plotfields_contours(els, xobs, yobs, 2, disp[:, 1], L"u_x")
    plotfields_contours(els, xobs, yobs, 3, disp[:, 2], L"u_y")
    plotfields_contours(els, xobs, yobs, 4, stress[:, 1], L"\sigma_{xx}")
    plotfields_contours(els, xobs, yobs, 5, stress[:, 2], L"\sigma_{yy}")
    plotfields_contours(els, xobs, yobs, 6, stress[:, 3], L"\sigma_{xy}")
    PyPlot.suptitle(suptitle, fontsize=20)
    show()
    return nothing
end

export rycroftcmap
"""
    rycroftcmap()

Chris Rycroft designed and shared this colormap.
"""
function rycroftcmap()
    cmap = ColorScheme([Colors.RGB(255.0 / 255.0, 20.0 / 255.0, 40.0 / 255.0),
                        Colors.RGB(255.0 / 255.0, 47.0 / 255.0, 146.0 / 255.0),
                        Colors.RGB(255.0 / 255.0, 138.0 / 255.0, 216.0 / 255.0),
                        Colors.RGB(55.0 / 255.0, 145.0 / 255.0, 230.0 / 255.0),
                        Colors.RGB(150.0 / 255.0, 230.0 / 255.0, 80.0 / 255.0),
                        Colors.RGB(255.0 / 255.0, 251.0 / 255.0, 0.0 / 255.0),
                        Colors.RGB(255.0 / 255.0, 255.0 / 255.0, 255.0 / 255.0)],
                        "rycroft", "psychedelic")
    cmap = ColorMap(ColorScheme([get(cmap, i) for i in 0.0:0.001:1.0]).colors)
    return cmap
end

export PUTC
"""
    PUTC(fun2dispstress, els, obsidx, srcidx, mu, nu)

Convenience function to calcualtion partial derivates of displacement,
and tractions for constant slip/displacement elements.  This
simply swaps the order of arguments (srcidx, obsidx) -> (obsidx, srcidx)
and does not return the stress partials.
"""
function PUTC(fun2dispstress, els, obsidx, srcidx, mu, nu)
    T, _, H = partialsconstdispstress(fun2dispstress, els, srcidx, obsidx, mu, nu)
    return T, H
end

export PUTQ
"""
    PUTQ(fun2dispstress, els, obsidx, srcidx, mu, nu)

Conveinece function to calcualtion partial derivates of displacement,
and tractions for constant slip/displacement elements.  This
simply swaps the order of arguments (srcidx, obsidx) -> (obsidx, srcidx)
and does not return the stress partials.
"""
function PUTQ(fun2dispstress, els, obsidx, srcidx, mu, nu)
    T, _, H = partialsquaddispstress(fun2dispstress, els, srcidx, obsidx, mu, nu)
    return T, H
end

export partialsconstdispstress
"""
    partialsconstdispstress(fun2dispstress, els, srcidx, obsidx, mu, nu)

Partial derivates of displacement, stress and tractions for constant
slip/displacement elements.

Note: Rename to parUSTC?
"""
function partialsconstdispstress(fun2dispstress, els, srcidx, obsidx, mu, nu)
    nobs, nsrc = length(obsidx), length(srcidx)
    partialsdisp, partialsstress, partialstrac = zeros(2 * nobs, 2 * nsrc), zeros(3 * nobs, 2 * nsrc), zeros(2 * nobs, 2 * nsrc)
    _partialsdisp, _partialsstress, _partialstrac = zeros(2, 2), zeros(3, 2), zeros(2, 2)
    for isrc in 1:nsrc
        for iobs in 1:nobs
            _partialsdisp[:, 1], _partialsstress[:, 1] = constdispstress(fun2dispstress, els.xcenter[obsidx[iobs]], els.ycenter[obsidx[iobs]], els, srcidx[isrc], 1, 0, mu, nu)
            _partialsdisp[:, 2], _partialsstress[:, 2] = constdispstress(fun2dispstress, els.xcenter[obsidx[iobs]], els.ycenter[obsidx[iobs]], els, srcidx[isrc], 0, 1, mu, nu)
            _partialstrac[:, 1] = stress2trac(_partialsstress[:, 1], [els.xnormal[obsidx[iobs]] ; els.ynormal[obsidx[iobs]]])
            _partialstrac[:, 2] = stress2trac(_partialsstress[:, 2], [els.xnormal[obsidx[iobs]] ; els.ynormal[obsidx[iobs]]])
            partialsdisp[2 * (iobs - 1) + 1:2 * (iobs - 1) + 2, 2 * (isrc - 1) + 1:2 * (isrc - 1) + 2] = _partialsdisp
            partialsstress[3 * (iobs - 1) + 1:3 * (iobs - 1) + 3, 2 * (isrc - 1) + 1:2 * (isrc - 1) + 2] = _partialsstress
            partialstrac[2 * (iobs - 1) + 1:2 * (iobs - 1) + 2, 2 * (isrc - 1) + 1:2 * (isrc - 1) + 2] = _partialstrac
        end
    end
    return partialsdisp, partialsstress, partialstrac
end

export quaddispstressinterleaved
"""
    quaddispstressinterleaved(fun2dispstress, xobs, yobs, els, idx, para, perp, mu, nu)

This is small helper function to return interleaved displacements
and stresses rather than the stacked columns.  For the non-coincident case.
"""
function quaddispstressinterleaved(fun2dispstress, xobs, yobs, els, idx, para, perp, mu, nu)
    dispstacked, stressstacked = quaddispstress(fun2dispstress, xobs, yobs, els, idx, para, perp, mu, nu)
    dispinterleaved = zeros(2 * length(xobs))
    dispinterleaved[1:2:end] = dispstacked[:, 1]
    dispinterleaved[2:2:end] = dispstacked[:, 2]
    stressinterleaved = zeros(3 * length(xobs))
    stressinterleaved[1:3:end] = stressstacked[:, 1]
    stressinterleaved[2:3:end] = stressstacked[:, 2]
    stressinterleaved[3:3:end] = stressstacked[:, 3]
    return dispinterleaved, stressinterleaved
end

export quaddispstresscoincidentinterleaved
"""
    quaddispstresscoincidentinterleaved(fun2dispstress, xobs, yobs, els, idx, para, perp, mu, nu)

This is small helper function to return interleaved displacements
and stresses rather than the stacked columns.  For the coincident case.
"""
function quaddispstresscoincidentinterleaved(fun2dispstress, xobs, yobs, els, idx, para, perp, mu, nu)
    dispstacked, stressstacked = quaddispstresscoincident(fun2dispstress, xobs, yobs, els, idx, para, perp, mu, nu)
    dispinterleaved = zeros(2 * length(xobs))
    dispinterleaved[1:2:end] = dispstacked[:, 1]
    dispinterleaved[2:2:end] = dispstacked[:, 2]
    stressinterleaved = zeros(3 * length(xobs))
    stressinterleaved[1:3:end] = stressstacked[:, 1]
    stressinterleaved[2:3:end] = stressstacked[:, 2]
    stressinterleaved[3:3:end] = stressstacked[:, 3]
    return dispinterleaved, stressinterleaved
end

export pointonline
"""
    pointonline(x1, y1, x2, y2, x3)

Given two endpointts of a line and an x value somewhere in between
calculate and return the corresponding y value
"""
function pointonline(x1, y1, x2, y2, x3)
    y3 = ((y2 - y1) * (x3 - x2)) / (x2 - x1) + y2
    return y3
end

export partialsquaddispstress
"""
    partialsquaddispstress(fun2dispstress, els, srcidx, obsidx, mu, nu)

Calculate partials derivaties for quadratic elements.

Each 6x6 displacement submatrix is shaped as:
       sx(ϕ1) sy(ϕ1) sx(ϕ2) sy(ϕ2) sx(ϕ3) sy(ϕ3)
ux(ϕ1)
uy(ϕ1)
ux(ϕ2)
uy(ϕ2)
ux(ϕ3)
uy(ϕ3)

Note: Rename to: parUSTQ?
"""
function partialsquaddispstress(fun2dispstress, els, srcidx, obsidx, mu, nu)
    nobs = length(obsidx)
    nsrc = length(srcidx)
    partialsdisp = zeros(6 * nobs, 6 * nsrc)
    partialsstress = zeros(9 * nobs, 6 * nsrc)
    partialstrac = zeros(6 * nobs, 6 * nsrc)
    _partialsdisp = zeros(6, 6)
    _partialsstress = zeros(9, 6)
    _partialstrac = zeros(6, 6)

    # Far-field interactions for all except main (block) diagonal
    for isrc in 1:nsrc
        for iobs in 1:nobs
            if srcidx[isrc] == obsidx[iobs]
                _partialsdisp[:, 1], _partialsstress[:, 1] = quaddispstresscoincidentinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [1 0 0], [0 0 0], mu, nu)
                _partialsdisp[:, 2], _partialsstress[:, 2] = quaddispstresscoincidentinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [0 0 0], [1 0 0], mu, nu)
                _partialsdisp[:, 3], _partialsstress[:, 3] = quaddispstresscoincidentinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [0 1 0], [0 0 0], mu, nu)
                _partialsdisp[:, 4], _partialsstress[:, 4] = quaddispstresscoincidentinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [0 0 0], [0 1 0], mu, nu)
                _partialsdisp[:, 5], _partialsstress[:, 5] = quaddispstresscoincidentinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [0 0 1], [0 0 0], mu, nu)
                _partialsdisp[:, 6], _partialsstress[:, 6] = quaddispstresscoincidentinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [0 0 0], [0 0 1], mu, nu)
            else
                _partialsdisp[:, 1], _partialsstress[:, 1] = quaddispstressinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [1 0 0], [0 0 0], mu, nu)
                _partialsdisp[:, 2], _partialsstress[:, 2] = quaddispstressinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [0 0 0], [1 0 0], mu, nu)
                _partialsdisp[:, 3], _partialsstress[:, 3] = quaddispstressinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [0 1 0], [0 0 0], mu, nu)
                _partialsdisp[:, 4], _partialsstress[:, 4] = quaddispstressinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [0 0 0], [0 1 0], mu, nu)
                _partialsdisp[:, 5], _partialsstress[:, 5] = quaddispstressinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [0 0 1], [0 0 0], mu, nu)
                _partialsdisp[:, 6], _partialsstress[:, 6] = quaddispstressinterleaved(fun2dispstress, els.xnodes[obsidx[iobs], :], els.ynodes[obsidx[iobs], :], els, srcidx[isrc], [0 0 0], [0 0 1], mu, nu)
            end
            for i in 1:6
                _partialstrac[1:2, i] = stress2trac(_partialsstress[1:3, i], [els.xnormal[obsidx[iobs]] ; els.ynormal[obsidx[iobs]]])
                _partialstrac[3:4, i] = stress2trac(_partialsstress[4:6, i], [els.xnormal[obsidx[iobs]] ; els.ynormal[obsidx[iobs]]])
                _partialstrac[5:6, i] = stress2trac(_partialsstress[7:9, i], [els.xnormal[obsidx[iobs]] ; els.ynormal[obsidx[iobs]]])
            end
            partialsdisp[6 * (iobs - 1) + 1:6 * (iobs - 1) + 6, 6 * (isrc - 1) + 1:6 * (isrc - 1) + 6] = _partialsdisp
            partialsstress[9 * (iobs - 1) + 1:9 * (iobs - 1) + 9, 6 * (isrc - 1) + 1:6 * (isrc - 1) + 6] = _partialsstress
            partialstrac[6 * (iobs - 1) + 1:6 * (iobs - 1) + 6, 6 * (isrc - 1) + 1:6 * (isrc - 1) + 6] = _partialstrac
        end
    end
    return partialsdisp, partialsstress, partialstrac
end

"""
    slip2coef(xnodes, slip, a)

Given a set of 3 displacements/slips at quadratic nodes convert
to a set of three phi values neccesary to comput the three quadratic
functions.
"""
function slip2coef(xnodes, slip, a)
    mat = zeros(3, 3)
    mat[:, 1] = (xnodes ./ a) .* (9 .* (xnodes ./ a) ./ 8 .- 3 ./ 4)
    mat[:, 2] = (1 .- 3 .* (xnodes ./ a) ./ 2) .* (1 .+ 3 .* (xnodes ./ a) ./ 2)
    mat[:, 3] = (xnodes ./ a) .* (9 .* (xnodes ./ a) ./ 8 .+ 3 ./ 4)
    return inv(mat) * slip
end

"""
    coef2slip(xnodes, slip, a)

Given a set of 3 phi values at quadratic nodes convert
to a set of three displacement/slip values neccesary to
physically interpret the results.
"""
function coef2slip(xnodes, coef, a)
    mat = zeros(3, 3)
    mat[:, 1] = (xnodes ./ a) .* (9 .* (xnodes ./ a) ./ 8 .- 3 ./ 4)
    mat[:, 2] = (1 .- 3 .* (xnodes ./ a) ./ 2) .* (1 .+ 3 * (xnodes ./ a) ./ 2)
    mat[:, 3] = (xnodes ./ a) .* (9 .* (xnodes ./ a) ./ 8 .+ 3 ./ 4)
    return mat * coef
end

export getidx
"""
    getidx(label, els)

Return a vector with the indices of els.name that match label.
"""
function getidx(label, els)
    idx = findall(x->x == label, els.name)
    println("getidx found " * string(length(idx)) * " elements with label \"" * label * "\"")
    return idx
end

export getidxdict
"""
    getidxdict(els)

Return a dictionary of the arrays containing the indices of els
associated with each distinct name in els.  The keys to the
dictionary are each unique name in els.
"""
function getidxdict(els)
    idxdict = Dict()
    names = unique(els.name)
    for i in 1:length(names)
        if length(names[i]) > 0
            idxdict[names[i]] = getidx(names[i], els)
        end
    end
    return idxdict
end

export initpartials
"""
    initpartials(els)

Create a nested dictionary for storing partial derivatives.
Should probably change this to Hstar, Sstar, and Tstar notation.
"""
function initpartials(els)
    partials = Dict()
    partials["disp"] = Dict()
    partials["stress"] = Dict()
    partials["trac"] = Dict()
    fieldnames = collect(keys(partials))
    elnames = collect(keys(getidxdict(els)))
    for iname in 1:length(fieldnames)
        for i in 1:length(elnames)
            partials[fieldnames[iname]][elnames[i]] = Dict()
            for j in 1:length(elnames)
                partials[fieldnames[iname]][elnames[i]][elnames[j]] = []
            end
        end
    end
    return partials
end

export quadstack
"""
    quadstack(vec)

Utility function for "stacking" a flattened quad vector
Useful for converting results of BEM solve to arguments
that can be used with quaddispstress for forward models
"""
function quadstack(vec)
    stack = transpose(reshape(vec, 3, Int(length(vec) / 3)))
    return stack
end

export plotqdtimeseries
"""
    plotqdtimeseries(sol, stridesize, nels)

Utility function for plotting quasidynaic results. For constant slip elements.
Can run out of memory :(
"""
function plotqdtimeseries(sol, stridesize, nels)
    siay = 365.25 * 24 * 60 * 60
    t = [x / siay for x in sol.t]
    theta = zeros(length(t), nels)
    vx = zeros(length(t), nels)
    vy = zeros(length(t), nels)
    for i in 1:length(t)
        vx[i, :] = sol.u[i][1:stridesize:end]
        vy[i, :] = sol.u[i][2:stridesize:end]
        theta[i, :] = sol.u[i][3:stridesize:end]
    end

    figure(1, figsize = (15, 8))
    subplot(3, 2, 1)
    plot(t, vx, "-", linewidth = 0.5)
    yscale("log")
    ylabel(L"v_x")

    subplot(3, 2, 3)
    plot(t, vy, "-", linewidth = 0.5)
    yscale("log")
    ylabel(L"v_y")

    subplot(3, 2, 5)
    plot(t, theta, "-", linewidth = 0.5)
    yscale("log")
    xlabel("t (years)")
    ylabel(L"\theta")

    subplot(3, 2, 2)
    plot(1:1:length(t), vx, "-", linewidth = 0.5)
    yscale("log")
    ylabel(L"v_x")

    subplot(3, 2, 4)
    plot(1:1:length(t), vy, "-", linewidth = 0.5)
    yscale("log")
    ylabel(L"v_y")

    subplot(3, 2, 6)
    plot(1:1:length(t), theta, "-", linewidth = 0.5)
    yscale("log")
    xlabel("time step #")
    ylabel(L"\theta")

    figure(2, figsize = (15, 5))
    plotme = log10.(vx')
    contourf(plotme, 5, cmap = rycroftcmap())
    colorbar(fraction = 0.020, pad = 0.05, extend = "both", label = L"$\log_{10}v$ (m/s)")
    contour(plotme, 5, linewidths = 0.5, linestyles = "solid", colors = "k")
    xlabel("time step")
    ylabel("element index")
    show()
end

export plotqdtimeseriesquad
"""
    plotqdtimeseriesquad(sol, stridesize, nels)

Utility function for plotting quasidynaic results. For quadratic slip elements.
Can run out of memory :(
"""
function plotqdtimeseriesquad(sol, stridesize, nels)
    siay = 365.25 * 24 * 60 * 60
    t = [x / siay for x in sol.t]
    theta = zeros(length(t), 3 * nels)
    vx = zeros(length(t), 3 * nels)
    vy = zeros(length(t), 3 * nels)
    for i in 1:length(t)
        vx[i, :] = sol.u[i][1:stridesize:end]
        vy[i, :] = sol.u[i][2:stridesize:end]
        theta[i, :] = sol.u[i][3:stridesize:end]
    end

    figure(figsize = (15, 8))
    subplot(3, 2, 1)
    plot(t, vx, "-", linewidth = 0.5)
    yscale("log")
    ylabel(L"v_x")

    subplot(3, 2, 3)
    plot(t, vy, "-", linewidth = 0.5)
    yscale("log")
    ylabel(L"v_y")

    subplot(3, 2, 5)
    plot(t, theta, "-", linewidth = 0.5)
    yscale("log")
    xlabel("t (years)")
    ylabel(L"\theta")

    subplot(3, 2, 2)
    plot(1:1:length(t), vx, "-", linewidth = 0.5)
    yscale("log")
    ylabel(L"v_x")

    subplot(3, 2, 4)
    plot(1:1:length(t), vy, "-", linewidth = 0.5)
    yscale("log")
    ylabel(L"v_y")

    subplot(3, 2, 6)
    plot(1:1:length(t), theta, "-", linewidth = 0.5)
    yscale("log")
    xlabel("time step #")
    ylabel(L"\theta")

    figure(figsize = (15, 5))
    plotme = log10.(vx')
    contourf(plotme, 5, cmap = rycroftcmap())
    colorbar(fraction = 0.020, pad = 0.05, extend = "both", label = L"$\log_{10}v$ (m/s)")
    contour(plotme, 5, linewidths = 0.5, linestyles = "solid", colors = "k")
    xlabel("time step")
    ylabel("element index")
    show()
end

export thetaaginglaw
"""
    thetaaginglaw(v, theta, dc)

Aging law in rate and state friction
"""
function thetaaginglaw(v, theta, dc)
    return @. 1 - theta * v / dc
end

export thetasliplaw
"""
    thetasliplaw(v, theta, dc)

Slip evolution law in rate and state friction
"""
function thetasliplaw(v, theta, dc)
    return @. -v * theta / dc * log(v * theta / dc)
end

export gravityparticularfunctions
"""
    gravityparticularfunctions()

From Pape and Bannerjee 1987
"""
function gravityparticularfunctions(x, y, g, rho, lambda, mu)
    U = zeros(length(x), 2)
    S = zeros(length(x), 3)
    U[:, 1] = @. -lambda * rho * g / (4 * mu * (lambda + mu)) * x * y  # Pape and Banerjee (1987) equation (6a)
    U[:, 2] = @. (rho * g) / (8 * mu * (lambda + mu)) * (lambda * x^2 + (lambda + 2 * mu) * y^2) # Pape and Banerjee (1987) equation (6b)
    S[:, 1] .= 0 # Pape and Banerjee (1987) equation (6c)
    S[:, 2] = @. rho * g * y  # Pape and Banerjee (1987) equation (6d)
    S[:, 3] .= 0 # Pape and Banerjee (1987) equation (6e)
    return U, S
end

end


text:
abstract type AbstractNormalSpline{n,T,RK} end

Base.ndims(::AbstractNormalSpline{n,T,RK}) where {n,T,RK} = n
Base.eltype(::AbstractNormalSpline{n,T,RK}) where {n,T,RK} = T

@inline _get_kernel(spl::AbstractNormalSpline)    = spl._kernel
@inline _get_nodes(spl::AbstractNormalSpline)     = spl._nodes
@inline _get_values(spl::AbstractNormalSpline)    = spl._values
@inline _get_d_nodes(spl::AbstractNormalSpline)   = spl._d_nodes
@inline _get_d_dirs(spl::AbstractNormalSpline)    = spl._d_dirs
@inline _get_d_values(spl::AbstractNormalSpline)  = spl._d_values
@inline _get_mu(spl::AbstractNormalSpline)        = spl._mu
@inline _get_rhs(spl::AbstractNormalSpline)       = spl._rhs
@inline _get_gram(spl::AbstractNormalSpline)      = spl._gram
@inline _get_chol(spl::AbstractNormalSpline)      = spl._chol
@inline _get_cond(spl::AbstractNormalSpline)      = spl._cond
@inline _get_min_bound(spl::AbstractNormalSpline) = spl._min_bound
@inline _get_max_bound(spl::AbstractNormalSpline) = spl._max_bound
@inline _get_scale(spl::AbstractNormalSpline)     = spl._scale

@doc raw"
`struct NormalSpline{n, T <: Real, RK <: ReproducingKernel_0} <: AbstractNormalSpline{n,T,RK}`

Define a structure containing full information of a normal spline
# Fields
- `_kernel`: a reproducing kernel spline was built with
- `_nodes`: transformed function value nodes
- `_values`: function values at interpolation nodes
- `_d_nodes`: transformed function directional derivative nodes
- `_d_dirs`: normalized derivative directions
- `_d_values`: function directional derivative values
- `_mu`: spline coefficients
- `_rhs`: right-hand side of the problem `gram * mu = rhs`
- `_gram`: Gram matrix of the problem `gram * mu = rhs`
- `_chol`: Cholesky factorization of the Gram matrix
- `_cond`: estimation of the Gram matrix condition number
- `_min_bound`: minimal bounds of the original node locations area
- `_max_bound`: maximal bounds of the original node locations area
- `_scale`: factor of transforming the original node locations into unit hypercube
"
Base.@kwdef struct NormalSpline{n, T <: Real, RK <: ReproducingKernel_0} <: AbstractNormalSpline{n,T,RK}
    _kernel::RK
    _nodes::VecOfSVecs{n,T}
    _values::Vector{T}        = zeros(T, 0)
    _d_nodes::VecOfSVecs{n,T} = zeros(SVector{n,T}, 0)
    _d_dirs::VecOfSVecs{n,T}  = zeros(SVector{n,T}, 0)
    _d_values::Vector{T}      = zeros(T, 0)
    _mu::Vector{T}            = zeros(T, 0)
    _rhs::Vector{T}           = zeros(T, 0)
    _gram::Hermitian{T, Matrix{T}}
    _chol::Cholesky{T, Matrix{T}}
    _cond::T
    _min_bound::SVector{n,T}
    _max_bound::SVector{n,T}
    _scale::T
end

Base.@kwdef struct ElasticNormalSpline{n, T <: Real, RK <: ReproducingKernel_0} <: AbstractNormalSpline{n,T,RK}
    _kernel::RK
    _max_size::Int
    _num_nodes::Base.RefValue{Int}      = Ref(0)
    _num_d_nodes::Base.RefValue{Int}    = Ref(0)
    _nodes::VecOfSVecs{n,T}             = zeros(SVector{n,T}, _max_size)
    _values::Vector{T}                  = zeros(T, _max_size)
    _d_nodes::VecOfSVecs{n,T}           = zeros(SVector{n,T}, n * _max_size)
    _d_dirs::VecOfSVecs{n,T}            = zeros(SVector{n,T}, n * _max_size)
    _d_values::Vector{T}                = zeros(T, n * _max_size)
    _mu::Vector{T}                      = zeros(T, (n+1) * _max_size)
    _rhs::Vector{T}                     = zeros(T, (n+1) * _max_size)
    _gram::Matrix{T}                    = zeros(T, (n+1) * _max_size, (n+1) * _max_size)
    _chol::ElasticCholesky{T,Matrix{T}} = ElasticCholesky{T}((n+1) * _max_size)
    _filled_columns::Vector{Int}        = zeros(Int, (n+1) * _max_size)
    _min_bound::SVector{n,T}
    _max_bound::SVector{n,T}
    _scale::T
end
function ElasticNormalSpline(min_bound::SVector{n,T}, max_bound::SVector{n,T}, max_size::Int, kernel::RK) where {n, T, RK <: ReproducingKernel_0}
    @assert kernel.ε != 0
    scale = maximum(max_bound .- min_bound)
    ElasticNormalSpline{n,T,RK}(; _kernel = kernel, _max_size = max_size, _min_bound = min_bound, _max_bound = max_bound, _scale = scale)
end

@inline _get_nodes(spl::ElasticNormalSpline)           = uview(spl._nodes, 1 : spl._num_nodes[])
@inline _get_values(spl::ElasticNormalSpline)          = uview(spl._values, 1 : spl._num_nodes[])
@inline _get_d_nodes(spl::ElasticNormalSpline)         = uview(spl._d_nodes, 1 : spl._num_d_nodes[])
@inline _get_d_dirs(spl::ElasticNormalSpline)          = uview(spl._d_dirs, 1 : spl._num_d_nodes[])
@inline _get_d_values(spl::ElasticNormalSpline)        = uview(spl._d_values, 1 : spl._num_d_nodes[])
@inline _get_cond(spl::ElasticNormalSpline)            = _estimate_cond(_get_gram(spl), _get_chol(spl))
@inline _get_mu(spl::ElasticNormalSpline)              = (J = _get_filled_columns(spl); return uview(spl._mu, J))
@inline _get_rhs(spl::ElasticNormalSpline)             = (J = _get_filled_columns(spl); return uview(spl._rhs, J))
@inline _get_gram(spl::ElasticNormalSpline)            = (J = _get_filled_columns(spl); A = uview(spl._gram, J, J); return Hermitian(A, :U))
@inline _get_filled_columns(spl::ElasticNormalSpline)  = uview(spl._filled_columns, 1 : spl._num_nodes[] + spl._num_d_nodes[])
@inline _get_insertion_order(spl::ElasticNormalSpline) = uview(_get_chol(spl).colperms, 1 : _get_chol(spl).ncols[])

function insertat!(x::AbstractVector, i, v, len = length(x))
    last = v
    @inbounds for j in i:len
        x[j], last = last, x[j]
    end
    return x
end

function Base.empty!(spl::ElasticNormalSpline)
    spl._num_nodes[] = 0
    spl._num_d_nodes[] = 0
    empty!(spl._chol)
    return spl
end

function Base.insert!(
        spl::ElasticNormalSpline{n,T,RK},
        node::SVector{n,T},
        value::T,
    ) where {n, T, RK <: ReproducingKernel_0}

    n₁, n₂ = spl._num_nodes[], spl._num_d_nodes[]
    n₁max  = spl._max_size
    n₂max  = n * spl._max_size
    @assert n₁ < n₁max

    # Normalize and insert node (assumed to be with `min_bound` and `_max_bound`)
    curr_nodes   = _get_nodes(spl)
    curr_d_nodes = _get_d_nodes(spl)
    curr_d_dirs  = _get_d_dirs(spl)
    new_node     = _normalize(spl, node)
    new_value    = value
    @inbounds begin
        spl._nodes[n₁+1]  = new_node
        spl._values[n₁+1] = new_value
        spl._rhs[n₁+1]    = new_value
        insertat!(spl._filled_columns, n₁+1, n₁+1, n₁+n₂+1)
        spl._num_nodes[] += 1
    end

    # Insert column into position `n₁+1` of Gram matrix
    inds = _get_filled_columns(spl)
    if RK <: ReproducingKernel_1
        _gram!(uview(spl._gram, inds, inds), new_node, curr_nodes, curr_d_nodes, curr_d_dirs, _get_kernel(spl))
    else
        _gram!(uview(spl._gram, inds, inds), new_node, curr_nodes, _get_kernel(spl))
    end

    # Insert column `n₁+1` of Gram matrix into Cholesky factorization
    insert!(spl._chol, n₁+1, Hermitian(spl._gram))
    cholesky!(spl._chol, n₁+1)

    # Solve for spline coefficients
    inds = _get_insertion_order(spl)
    ldiv!(uview(spl._mu, inds), spl._chol, uview(spl._rhs, inds))

    return nothing
end

function Base.insert!(
        spl::ElasticNormalSpline{n,T,RK},
        d_node::SVector{n},
        d_dir::SVector{n},
        d_value::T,
    ) where {n, T, RK <: ReproducingKernel_1}

    n₁, n₂ = spl._num_nodes[], spl._num_d_nodes[]
    n₁max  = spl._max_size
    n₂max  = n * spl._max_size
    @assert n₂ < n₂max

    # Normalize and insert node (assumed to be with `min_bound` and `_max_bound`)
    curr_nodes   = _get_nodes(spl)
    curr_d_nodes = _get_d_nodes(spl)
    curr_d_dirs  = _get_d_dirs(spl)
    new_d_node   = _normalize(spl, d_node)
    new_d_dir    = d_dir / norm(d_dir)
    new_d_value  = _get_scale(spl) * d_value
    @inbounds begin
        spl._d_nodes[n₂+1]   = new_d_node
        spl._d_values[n₂+1]  = new_d_value
        spl._d_dirs[n₂+1]    = new_d_dir
        spl._rhs[n₁max+n₂+1] = new_d_value
        insertat!(spl._filled_columns, n₁+n₂+1, n₁max+n₂+1, n₁+n₂+1)
        spl._num_d_nodes[]  += 1
    end

    # Insert column into position `n₁max+n₂+1` of Gram matrix
    inds = _get_filled_columns(spl)
    _gram!(uview(spl._gram, inds, inds), new_d_node, new_d_dir, curr_nodes, curr_d_nodes, curr_d_dirs, _get_kernel(spl))

    # Insert column `n₁max+n₂+1` of Gram matrix into Cholesky factorization
    insert!(spl._chol, n₁max+n₂+1, Hermitian(spl._gram))
    cholesky!(spl._chol, n₁max+n₂+1)

    # Solve for spline coefficients
    inds = _get_insertion_order(spl)
    ldiv!(uview(spl._mu, inds), spl._chol, uview(spl._rhs, inds))

    return nothing
end

function Base.insert!(
        spl::ElasticNormalSpline{n,T,RK},
        nodes::AbstractVecOfSVecs{n,T},
        values::AbstractVector{T},
    ) where {n, T, RK <: ReproducingKernel_0}
    @assert length(nodes) == length(values)

    # Insert `n` regular nodes
    @inbounds for i in 1:length(nodes)
        insert!(spl, nodes[i], values[i])
    end

    return spl
end

function Base.insert!(
        spl::ElasticNormalSpline{n,T,RK},
        nodes::AbstractVecOfSVecs{n,T},
        values::AbstractVector{T},
        d_nodes::AbstractVecOfSVecs{n,T},
        d_dirs::AbstractVecOfSVecs{n,T},
        d_values::AbstractVector{T},
    ) where {n, T, RK <: ReproducingKernel_1}
    @assert length(nodes) == length(values)
    @assert length(d_nodes) == length(d_dirs) == length(d_values)

    # Insert `n` regular nodes
    @inbounds for i in 1:length(nodes)
        insert!(spl, nodes[i], values[i])
    end

    # Insert `n` derivative nodes
    @inbounds for i in 1:length(d_nodes)
        insert!(spl, d_nodes[i], d_dirs[i], d_values[i])
    end

    return spl
end

ext: jl


text:
# # Extending `@set` and `@lens`
# This code demonstrates how to extend the `@set` and `@lens` mechanism with custom
# lenses.
# As a demo, we want to implement `@mylens!` and `@myset!`, which work much like
# `@lens` and `@set`, but mutate objects instead of returning modified copies.

using Setfield
using Setfield: IndexLens, PropertyLens, ComposedLens

struct Lens!{L <:Lens} <: Lens
    pure::L
end

Setfield.get(o, l::Lens!) = Setfield.get(o, l.pure)
function Setfield.set(o, l::Lens!{<: ComposedLens}, val)
    o_inner = get(o, l.pure.outer)
    set(o_inner, Lens!(l.pure.inner), val)
end
function Setfield.set(o, l::Lens!{PropertyLens{prop}}, val) where {prop}
    setproperty!(o, prop, val)
    o
end
function Setfield.set(o, l::Lens!{<:IndexLens}, val) where {prop}
    o[l.pure.indices...] = val
    o
end

# Now this implements the kind of `lens` the new macros should use.
# Of course there are more variants like `Lens!(<:DynamicIndexLens)`, for which we might
# want to overload `set`, but lets ignore that. Instead we want to check, that everything works so far:

using Test
mutable struct M
    a
    b
end

o = M(1,2)
l = Lens!(@lens _.b)
set(o, l, 20)
@test o.b == 20

l = Lens!(@lens _.foo[1])
o = (foo=[1,2,3], bar=:bar)
set(o, l, 100)
@test o == (foo=[100,2,3], bar=:bar)

# Now we can implement the syntax macros

using Setfield: setmacro, lensmacro

macro myset!(ex)
    setmacro(Lens!, ex)
end

macro mylens!(ex)
    lensmacro(Lens!, ex)
end

o = M(1,2)
@myset! o.a = :hi
@myset! o.b += 98
@test o.a == :hi
@test o.b == 100

deep = [[[[1]]]]
@myset! deep[1][1][1][1] = 2
@test deep[1][1][1][1] === 2

l = @mylens! _.foo[1]
o = (foo=[1,2,3], bar=:bar)
set(o, l, 100)
@test o == (foo=[100,2,3], bar=:bar)

# Everything works, we can do arbitrary nesting and also use `+=` syntax etc.


text:
using Base.Test

using ucbature

const my_tol = 1e-2

const dim = 2
const nb_evals = 1e6
const xmin = fill(-1.0, dim)
const xmax = fill( 1.0, dim)
const nb_zones_per_dim = 5
const k_exploration = 0.01
const nb_init_evals_per_zone = 10

sphere(x) = vecdot(x, x) < 1.0 ? 1.0 : 0.0

res_mc = ucbature.mc(sphere, xmin, xmax, nb_evals)
@test isapprox(res_mc, pi, atol=my_tol)

res_str = ucbature.str(sphere, xmin, xmax, nb_evals, nb_zones_per_dim)
@test isapprox(res_str, pi, atol=my_tol)

res_ucb = ucbature.ucb(sphere, xmin, xmax, nb_evals, nb_zones_per_dim,
                       k_exploration, nb_init_evals_per_zone)
@test isapprox(res_ucb, pi, atol=my_tol)

max_stars_repo_name: florianLepretre/ucbature

max_issues_repo_name: florianLepretre/ucbature

max_forks_repo_name: florianLepretre/ucbature


text:
""" Computer algebra via monoidal categories.

In a conventional computer algebra system, algebraic expressions are represented
as *trees* whose leaves are variables or constants and whose internal nodes are
arithmetic operations or elementary or special functions. The idea here is to
represent expressions as morphisms in a monoidal category.
"""
module AlgebraicNets
export AlgebraicNetTheory, AlgebraicNet, Ob, Hom, dom, codom,
  id, compose, ⋅, ∘, otimes, ⊗, munit, braid, mcopy, delete, mmerge, create,
  linear, constant,
  compile, compile_expr, compile_expr_vector, compile_block, evaluate

using GeneralizedGenerated: mk_function
using Match
import StaticArrays

using Catlab
using Catlab.Meta: concat_expr
import Catlab.Syntax: show_latex, show_unicode
using Catlab.Theories: MonoidalCategoryWithBidiagonals, ObExpr, HomExpr
import Catlab.Theories: Ob, Hom, dom, codom,
  id, compose, ⋅, ∘, otimes, ⊗, munit, braid, mcopy, delete, mmerge, create
using Catlab.Programs
import Catlab.Programs: compile, compile_expr, compile_block, evaluate_hom
import Catlab.Programs.GenerateJuliaPrograms: genvar, genvars, to_function_expr,
  generator_expr, input_exprs

# Syntax
########

""" Theory of *algebraic networks*

TODO: Explain
"""
@signature AlgebraicNetTheory{Ob,Hom} <: MonoidalCategoryWithBidiagonals{Ob,Hom} begin
  linear(x::Any, A::Ob, B::Ob)::(A → B)
  constant(x::Any, A::Ob)::(munit() → A)
end

@syntax AlgebraicNet{ObExpr,HomExpr} AlgebraicNetTheory begin
  # FIXME: `compose` and `otimes` should delegate to wiring layer when possible.
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))
  otimes(A::Ob, B::Ob) = associate_unit(new(A,B), munit)
  otimes(f::Hom, g::Hom) = associate(new(f,g))

  mcopy(A::Ob) = mcopy(A, 2)
  mmerge(A::Ob) = mmerge(A, 2)
end

function mcopy(A::AlgebraicNet.Ob, n::Int)
  AlgebraicNet.Hom{:mcopy}([A, n], [A, otimes(fill(A, n))])
end
function mmerge(A::AlgebraicNet.Ob, n::Int)
  AlgebraicNet.Hom{:mmerge}([A, n], [otimes(fill(A, n)), A])
end

# Compilation
#############

""" Internal state for compilation of algebraic network into Julia function.
"""
mutable struct AlgebraicNetState <: CompileState
  nvars::Int
  constants::Dict{Symbol,Int}
  constants_sym::Symbol
  AlgebraicNetState(; nvars=0, constants=Dict{Symbol,Int}(), constants_sym=Symbol()) =
    new(nvars, constants, constants_sym)
end

""" Compile an algebraic network into a Julia function.

This method of "functorial compilation" generates simple imperative code with no
optimizations. Still, the code should be fast provided the original expression
is properly factored, with no duplicate computations.
"""
function compile(f::Union{AlgebraicNet.Hom,Block};
                 return_constants::Bool=false, vector::Bool=false, kw...)
  expr, constants = vector ? compile_expr_vector(f; kw...) : compile_expr(f; kw...)
  compiled = mk_function(Main, expr)
  return_constants ? (compiled, constants) : compiled
end

""" Compile an algebraic network into a Julia function expression.

The function signature is:
  - arguments = inputs (domain) of network
  - keyword arguments = symbolic constants (coefficients) of network, if any
"""
function compile_expr(f::AlgebraicNet.Hom; args::Vector{Symbol}=Symbol[])
  inputs = isempty(args) ? input_exprs(ndims(dom(f)), kind=:variables) : args
  block, constants = compile_block(f, inputs)
  function_expr = to_function_expr(block; kwargs=constants)
  (function_expr, constants)
end

""" Compile an algebraic network into a Julia function expression.

The function signature is:
  - first argument = input vector
  - second argument = constant (coefficients) vector

Unlike `compile_expr`, this method assumes the network has a single output.
"""
function compile_expr_vector(f::AlgebraicNet.Hom; inputs_sym::Symbol=:x,
                             constants_sym::Symbol=:c)
  # Compile algebraic network in block.
  inputs = input_exprs(ndims(dom(f)), prefix=inputs_sym, kind=:array)
  block, constants = compile_block(f, inputs, constants_sym=constants_sym)

  # Create call expression (function header).
  call_expr = Expr(:tuple, inputs_sym, constants_sym)

  # Create function body.
  @assert length(block.outputs) == 1
  return_expr = Expr(:return, block.outputs[1])
  body_expr = concat_expr(block.code, return_expr)

  (Expr(:function, call_expr, body_expr), constants)
end

""" Compile an algebraic network into a block of Julia code.
"""
function compile_block(f::AlgebraicNet.Hom, inputs::Vector; kw...)::Tuple{Block,Vector}
  state = AlgebraicNetState(; kw...)
  block = compile_block(f, inputs, state)
  constants = [ k for (k,v) in sort(collect(state.constants), by=x->x[2]) ]
  return (block, constants)
end

function compile_block(f::AlgebraicNet.Hom{:linear}, inputs::Vector,
                       state::AlgebraicNetState)::Block
  nin, nout = ndims(dom(f)), ndims(codom(f))
  outputs = genvars(state, nout)
  @assert length(inputs) == nin

  value = first(f)
  if value isa Symbol
    value = genconst(state, value)
  elseif value isa AbstractMatrix
    value = Expr(:call,
      Expr(:curly, GlobalRef(StaticArrays, :SMatrix), size(value)...),
      value...)
  end
  lhs = nout == 1 ? outputs[1] : Expr(:tuple, outputs...)
  rhs = if nin == 0
    0
  elseif nin == 1
    Expr(:call, :(*), value, inputs[1])
  else
    Expr(:call, :(*), value, Expr(:vect, inputs...))
  end
  Block(Expr(:(=), lhs, rhs), inputs, outputs)
end

function compile_block(f::AlgebraicNet.Hom{:constant}, inputs::Vector,
                       state::AlgebraicNetState)::Block
  @assert isempty(inputs)
  outputs = genvars(state, ndims(codom(f)))
  value = first(f)
  value = isa(value, Symbol) ? genconst(state, value) : value
  lhs = length(outputs) == 1 ? outputs[1] : Expr(:tuple, outputs...)
  Block(Expr(:(=), lhs, value), inputs, outputs)
end

# FIXME: Wiring layers are gone. Migrate to another data structure?
#=
function compile_block(f::AlgebraicNet.Hom{:wiring}, inputs::Vector,
                       state::CompileState)::Block
  nout = ndims(codom(f))
  outputs = genvars(state, nout)
  terms = [ [] for i in 1:nout ]
  for (src, tgts) in first(f).wires
    x = inputs[src]
    for (tgt, c) in tgts
      push!(terms[tgt], c == 1 ? x : Expr(:call, :(.*), c, x))
    end
  end
  code = multiple_assign_expr(outputs, map(sum_expr, terms))
  Block(code, inputs, outputs)
end
=#

function compile_block(f::AlgebraicNet.Hom{:mmerge}, inputs::Vector,
                       state::CompileState)::Block
  out = genvar(state)
  code = Expr(:(=), out, sum_expr(inputs))
  Block(code, inputs, [out])
end

function compile_block(f::AlgebraicNet.Hom{:create}, inputs::Vector,
                       state::CompileState)::Block
  nout = ndims(codom(f))
  outputs = genvars(state, nout)
  code = multiple_assign_expr(outputs, zeros(nout))
  Block(code, inputs, outputs)
end

""" Generate Julia expression for morphism generator in algebraic net.
"""
function generator_expr(f::AlgebraicNet.Hom{:generator}, inputs::Vector,
                        state::AlgebraicNetState)
  value = first(f)
  if isempty(inputs)
    # If nullary, a constant, possibly symbolic.
    isa(value, Symbol) ? genconst(state, value) : value
  else
    # Otherwise, a broadcasting function call.
    Expr(:(.), value, Expr(:tuple, inputs...))
  end
end

""" Generate Julia expression for single or multiple assignment.
"""
function multiple_assign_expr(lhs::Vector, rhs::Vector)::Expr
  @assert length(lhs) == length(rhs)
  if length(lhs) == 1
    Expr(:(=), lhs[1], rhs[1])
  else
    Expr(:(=), Expr(:tuple, lhs...), Expr(:tuple, rhs...))
  end
end

""" Generate Julia expression for sum of zero, one, or more terms.
"""
function sum_expr(T::Type, terms::Vector)
  if length(terms) == 0
    zero(T)
  elseif length(terms) == 1
    terms[1]
  else
    Expr(:call, GlobalRef(Base, :broadcast), :(+), terms...)
  end
end
sum_expr(terms::Vector) = sum_expr(Float64, terms)

""" Generate a constant (symbol or expression).

See also `gensym` and `genvar`.
"""
function genconst(state::AlgebraicNetState, name::Symbol)
  i = get!(state.constants, name, length(state.constants)+1)
  state.constants_sym == Symbol() ? name : :($(state.constants_sym)[$i])
end

# Evaluation
############

function evaluate_hom(f::AlgebraicNet.Hom{:generator}, xs::Vector; kw...)
  value = first(f)
  y = if value isa Symbol && head(dom(f)) != :munit
    # Broadcasting function call.
    getfield(Main, value).(xs...)
  else
    value
  end
  y isa Tuple ? collect(y) : [y]
end

function evaluate_hom(f::AlgebraicNet.Hom{:wiring}, xs::Vector; kw...)
  # XXX: We can't properly preallocate the y's because we don't their dims.
  ys = repeat(Any[0.0], ndims(codom(f)))
  for (src, tgts) in first(f).wires
    for (tgt, c) in tgts
      y = c * xs[src]
      ys[tgt] == 0 ? (ys[tgt] = y) : (ys[tgt] += y)
    end
  end
  ys
end

evaluate_hom(f::AlgebraicNet.Hom{:mmerge}, xs::Vector; kw...) = [.+(xs...)]
evaluate_hom(f::AlgebraicNet.Hom{:create}, xs::Vector; kw...) = zeros(ndims(codom(f)))
evaluate_hom(f::AlgebraicNet.Hom{:linear}, xs::Vector; kw...) = first(f) * xs
evaluate_hom(f::AlgebraicNet.Hom{:constant}, xs::Vector; kw...) = first(f)

# Display
#########

""" Denote composition by a semicolon ala the computer scientists.

In this context, `⋅` is too easily confused for multiplication, ` ` (space) is
too implicit, and `∘` has a right-to-left connotation.
"""
function show_latex(io::IO, expr::AlgebraicNet.Hom{:compose}; paren::Bool=false, kw...)
  Syntax.show_latex_infix(io, expr, ";"; paren=paren, kw...)
end
function show_unicode(io::IO, expr::AlgebraicNet.Hom{:compose}; kw...)
  Syntax.show_unicode_infix(io, expr, "; "; kw...)
end

function show_latex(io::IO, expr::AlgebraicNet.Hom{:linear}; kw...)
  print(io, "\\mathop{\\mathrm{linear}}\\left[$(first(expr))\\right]")
end
function show_unicode(io::IO, expr::AlgebraicNet.Hom{:linear}; kw...)
  print(io, "linear{$(first(expr))}")
end
function show_latex(io::IO, expr::AlgebraicNet.Hom{:constant}; kw...)
  print(io, first(expr))
end
function show_unicode(io::IO, expr::AlgebraicNet.Hom{:constant}; kw...)
  print(io, first(expr))
end

end


text:
module PQBaseCamp

# include my codes -
include("Include.jl")

# functions/types that we are going to return -
export Δ
export 𝒟
export 𝒫
export β
export covariance
export sample

# export models -
export AbstractBaseCampComputation
export LogReturnComputationModel
export LinearReturnComputationModel

end # module


text:
immutable NoncentralT <: ContinuousUnivariateDistribution
    df::Float64
    ncp::Float64
    function NonCentralT(d::Real, nc::Real)
    	if d >= 0.0 && nc >= 0.0
    		new(float64(d), float64(nc))
    	else
    		error("df and ncp must be non-negative")
    	end
    end
end

@_jl_dist_2p NoncentralT nt

insupport(d::NoncentralT, x::Number) = isreal(x) && isfinite(x)


text:
using Documenter, AMDGPUnative

makedocs(
    sitename="AMDGPUnative.jl",
    pages = [
        "Home" => "index.md",
        "Quick Start" => "quickstart.md",
        "Global Variables" => "globals.md",
        "API Reference" => "api.md"
    ]
)


text:
# ## `items/image.jl`
#
# ### Item
#
# We first define the [`Image`](#) item. Since we need to keep
# track of the spatial bounds for projective transformations
# we add them as a field. By default, they will simply
# correspond to the image axes.

"""
    Image(image[, bounds])

Item representing an N-dimensional image with element type T.

## Examples

{cell=image}
```julia
using DataAugmentation, Images

imagedata = rand(RGB, 100, 100)
item = Image(imagedata)
showitem(item)
```

If `T` is not a color, the image will be interpreted as grayscale:

{cell=image}
```julia
imagedata = rand(Float32, 100, 100)
item = Image(imagedata)
showitem(item)
```

"""
struct Image{N,T,B} <: AbstractArrayItem{N,T}
    data::AbstractArray{T,N}
    bounds::AbstractArray{<:SVector{N,B},N}
end

Image(data) = Image(data, size(data))

function Image(data::AbstractArray{T,N}, sz::NTuple{N,Int}) where {T,N}
    bounds = makebounds(sz)
    return Image(data, bounds)
end

Base.show(io::IO, item::Image{N,T}) where {N,T} =
    print(io, "Image{$N, $T}() with size $(size(itemdata(item)))")

function showitem(image::Image{2, <:Colorant})
    return image.data
end

function showitem(image::Image{2, <:AbstractFloat})
    return colorview(Gray, image.data)
end

# ### Projective transformations
#
# To support projective transformations, we need to implement [`getbounds`](#)
# and [`project`](#).

getbounds(image::Image) = image.bounds

# For the projection, we use `warp` from
# [ImageTransformations.jl](https://github.com/JuliaImages/ImageTransformations.jl).
# We have to pass the inverse of the projection `P` as it uses backward
# mode warping.

function project(P, image::Image{N, T}, indices) where {N, T}
    ## Transform the bounds along with the image
    bounds_ = P.(getbounds(image))
    data_ = warp(itemdata(image), inv(P), indices, zero(T))
    return Image(data_, makebounds(indices))
end

# The inplace version `project!` is quite similar. Note `indices` are not needed
# as they are implicitly given by the buffer.

function project!(bufimage::Image, P, image::Image{N, T}, indices) where {N, T}
    warp!(
        itemdata(bufimage),
        box_extrapolation(itemdata(image), zero(T)),
        inv(P),
    )
    return Image(itemdata(bufimage), P.(getbounds(image)))
end


text:
"""
    VerstappenAnisotropicMinimumDissipation{T} <: AbstractAnisotropicMinimumDissipation{T}

Use the anisotropic minimum dissipation large eddy simulation model proposed by Verstappen
(2018) and described by Vreugdenhil & Taylor (2018).
"""
struct VerstappenAnisotropicMinimumDissipation{T} <: AbstractAnisotropicMinimumDissipation{T}
     C :: T
    Cb :: T
     ν :: T
     κ :: T
end

"""
    VerstappenAnisotropicMinimumDissipation(T=Float64; C = 1/12, Cb = 0.0, ν = ν₀, κ=κ₀)

Use a `VerstappenAnisotropicMinimumDissipation` turbulence closure with parameters of type
`T`.

Keyword arguments
=================
    - `C::T`: Poincaré constant
    - `Cb::T`: Buoyancy modification multiplier (`Cb = 0` turns it off, `Cb = 1` turns it on)
    - `ν::T`: 'molecular' background viscosity for momentum
    - `κ::T`: 'molecular' background diffusivity for tracers

By default, `C` = 1/12 which is appropriate for a finite-volume method employing a
second-order advection scheme, `Cb` = 0 which terms off the buoyancy modification term,
and molecular values are used for `ν` and `κ`.

References
==========
Vreugdenhil C., and Taylor J. (2018), "Large-eddy simulations of stratified plane Couette
    flow using the anisotropic minimum-dissipation model", Physics of Fluids 30, 085104.
Verstappen, R. (2018), "How much eddy dissipation is needed to counterbalance the nonlinear
    production of small, unresolved scales in a large-eddy simulation of turbulence?",
    Computers & Fluids 176, pp. 276-284.
"""
function VerstappenAnisotropicMinimumDissipation(T=Float64;
     C = 1/12,
    Cb = 0.0,
     ν = ν₀,
     κ = κ₀,
    )
    return VerstappenAnisotropicMinimumDissipation{T}(C, Cb, ν, κ)
end

const VAMD = VerstappenAnisotropicMinimumDissipation

function TurbulentDiffusivities(arch::AbstractArchitecture, grid::AbstractGrid, ::VAMD)
     νₑ = CellField(arch, grid)
    κTₑ = CellField(arch, grid)
    κSₑ = CellField(arch, grid)
    return (νₑ=νₑ, κₑ=(T=κTₑ, S=κSₑ))
end

@inline function ν_ccc(i, j, k, grid::AbstractGrid{FT}, closure::VAMD, c, buoyancy, U, Φ) where FT
    ijk = (i, j, k, grid)
    q = norm_tr_∇u_ccc(ijk..., U.u, U.v, U.w)

    if q == 0 # SGS viscosity is zero when strain is 0
        νˢᵍˢ = zero(FT)
    else
        r = norm_uᵢₐ_uⱼₐ_Σᵢⱼ_ccc(ijk..., closure, U.u, U.v, U.w)
        ζ = norm_wᵢ_bᵢ_ccc(ijk..., closure, buoyancy, U.w, Φ) / Δᶠz_ccc(ijk...)
        δ² = 3 / (1 / Δᶠx_ccc(ijk...)^2 + 1 / Δᶠy_ccc(ijk...)^2 + 1 / Δᶠz_ccc(ijk...)^2)
        νˢᵍˢ = - closure.C * δ² * (r - closure.Cb * ζ) / q
    end

    return max(zero(FT), νˢᵍˢ) + closure.ν
end

@inline function κ_ccc(i, j, k, grid::AbstractGrid{FT}, closure::VAMD, c, buoyancy, U, Φ) where FT
    ijk = (i, j, k, grid)
    σ =  norm_θᵢ²_ccc(i, j, k, grid, c) # Tracer variance

    if σ == 0
        κˢᵍˢ = zero(FT)
    else
        ϑ =  norm_uᵢⱼ_cⱼ_cᵢ_ccc(ijk..., closure, U.u, U.v, U.w, c)
        δ² = 3 / (1 / Δᶠx_ccc(ijk...)^2 + 1 / Δᶠy_ccc(ijk...)^2 + 1 / Δᶠz_ccc(ijk...)^2)
        κˢᵍˢ = - closure.C * δ² * ϑ / σ
    end

    return max(zero(FT), κˢᵍˢ) + closure.κ
end

#####
##### The *** 30 terms *** of AMD
#####

@inline function norm_uᵢₐ_uⱼₐ_Σᵢⱼ_ccc(i, j, k, grid, closure, u, v, w)
    ijk = (i, j, k, grid)
    uvw = (u, v, w)
    ijkuvw = (i, j, k, grid, u, v, w)

    uᵢ₁_uⱼ₁_Σ₁ⱼ = (
         norm_Σ₁₁(ijkuvw...) * norm_∂x_u(ijk..., u)^2
      +  norm_Σ₂₂(ijkuvw...) * ▶xy_cca(ijk..., norm_∂x_v², uvw...)
      +  norm_Σ₃₃(ijkuvw...) * ▶xz_cac(ijk..., norm_∂x_w², uvw...)

      +  2 * norm_∂x_u(ijkuvw...) * ▶xy_cca(ijk..., norm_∂x_v_Σ₁₂, uvw...)
      +  2 * norm_∂x_u(ijkuvw...) * ▶xz_cac(ijk..., norm_∂x_w_Σ₁₃, uvw...)
      +  2 * ▶xy_cca(ijk..., norm_∂x_v, uvw...) * ▶xz_cac(ijk..., norm_∂x_w, uvw...)
           * ▶yz_acc(ijk..., norm_Σ₂₃, uvw...)
    )

    uᵢ₂_uⱼ₂_Σ₂ⱼ = (
      + norm_Σ₁₁(ijkuvw...) * ▶xy_cca(ijk..., norm_∂y_u², uvw...)
      + norm_Σ₂₂(ijkuvw...) * norm_∂y_v(ijk..., v)^2
      + norm_Σ₃₃(ijkuvw...) * ▶yz_acc(ijk..., norm_∂y_w², uvw...)

      +  2 * norm_∂y_v(ijkuvw...) * ▶xy_cca(ijk..., norm_∂y_u_Σ₁₂, uvw...)
      +  2 * ▶xy_cca(ijk..., norm_∂y_u, uvw...) * ▶yz_acc(ijk..., norm_∂y_w, uvw...)
           * ▶xz_cac(ijk..., norm_Σ₁₃, uvw...)
      +  2 * norm_∂y_v(ijkuvw...) * ▶yz_acc(ijk..., norm_∂y_w_Σ₂₃, uvw...)
    )

    uᵢ₃_uⱼ₃_Σ₃ⱼ = (
      + norm_Σ₁₁(ijkuvw...) * ▶xz_cac(ijk..., norm_∂z_u², uvw...)
      + norm_Σ₂₂(ijkuvw...) * ▶yz_acc(ijk..., norm_∂z_v², uvw...)
      + norm_Σ₃₃(ijkuvw...) * norm_∂z_w(ijk..., w)^2

      +  2 * ▶xz_cac(ijk..., norm_∂z_u, uvw...) * ▶yz_acc(ijk..., norm_∂z_v, uvw...)
           * ▶xy_cca(ijk..., norm_Σ₁₂, uvw...)
      +  2 * norm_∂z_w(ijkuvw...) * ▶xz_cac(ijk..., norm_∂z_u_Σ₁₃, uvw...)
      +  2 * norm_∂z_w(ijkuvw...) * ▶yz_acc(ijk..., norm_∂z_v_Σ₂₃, uvw...)
    )

    return uᵢ₁_uⱼ₁_Σ₁ⱼ + uᵢ₂_uⱼ₂_Σ₂ⱼ + uᵢ₃_uⱼ₃_Σ₃ⱼ
end

#####
##### trace(∇u) = uᵢⱼ uᵢⱼ
#####

@inline function norm_tr_∇u_ccc(i, j, k, grid, uvw...)
    ijk = (i, j, k, grid)

    return (
        # ccc
        norm_∂x_u²(ijk..., uvw...)
      + norm_∂y_v²(ijk..., uvw...)
      + norm_∂z_w²(ijk..., uvw...)

        # ffc
      + ▶xy_cca(ijk..., norm_∂x_v², uvw...)
      + ▶xy_cca(ijk..., norm_∂y_u², uvw...)

        # fcf
      + ▶xz_cac(ijk..., norm_∂x_w², uvw...)
      + ▶xz_cac(ijk..., norm_∂z_u², uvw...)

        # cff
      + ▶yz_acc(ijk..., norm_∂y_w², uvw...)
      + ▶yz_acc(ijk..., norm_∂z_v², uvw...)
    )
end

@inline function norm_wᵢ_bᵢ_ccc(i, j, k, grid, closure, buoyancy, w, C)
    ijk = (i, j, k, grid)

    wx_bx = (▶xz_cac(ijk..., norm_∂x_w, w)
             * Δᶠx_ccc(ijk...) * ▶x_caa(ijk..., ∂x_faa, buoyancy_perturbation, buoyancy, C))

    wy_by = (▶yz_acc(ijk..., norm_∂y_w, w)
             * Δᶠy_ccc(ijk...) * ▶y_aca(ijk..., ∂y_afa, buoyancy_perturbation, buoyancy, C))

    wz_bz = (norm_∂z_w(ijk..., w)
             * Δᶠz_ccc(ijk...) * ▶z_aac(ijk..., ∂z_aaf, buoyancy_perturbation, buoyancy, C))

    return wx_bx + wy_by + wz_bz
end

@inline function norm_uᵢⱼ_cⱼ_cᵢ_ccc(i, j, k, grid, closure, u, v, w, c)
    ijk = (i, j, k, grid)

    cx_ux = (
                   norm_∂x_u(ijk..., u) * ▶x_caa(ijk..., norm_∂x_c², c)
        + ▶xy_cca(ijk..., norm_∂x_v, v) * ▶x_caa(ijk..., norm_∂x_c, c) * ▶y_aca(ijk..., norm_∂y_c, c)
        + ▶xz_cac(ijk..., norm_∂x_w, w) * ▶x_caa(ijk..., norm_∂x_c, c) * ▶z_aac(ijk..., norm_∂z_c, c)
    )

    cy_uy = (
          ▶xy_cca(ijk..., norm_∂y_u, u) * ▶y_aca(ijk..., norm_∂y_c, c) * ▶x_caa(ijk..., norm_∂x_c, c)
        +          norm_∂y_v(ijk..., v) * ▶y_aca(ijk..., norm_∂y_c², c)
        + ▶xz_cac(ijk..., norm_∂y_w, w) * ▶y_aca(ijk..., norm_∂y_c, c) * ▶z_aac(ijk..., norm_∂z_c, c)
    )

    cz_uz = (
          ▶xz_cac(ijk..., norm_∂z_u, u) * ▶z_aac(ijk..., norm_∂z_c, c) * ▶x_caa(ijk..., norm_∂x_c, c)
        + ▶yz_acc(ijk..., norm_∂z_v, v) * ▶z_aac(ijk..., norm_∂z_c, c) * ▶y_aca(ijk..., norm_∂y_c, c)
        +          norm_∂z_w(ijk..., w) * ▶z_aac(ijk..., norm_∂z_c², c)
    )

    return cx_ux + cy_uy + cz_uz
end

@inline norm_θᵢ²_ccc(i, j, k, grid, c) = (
      ▶x_caa(i, j, k, grid, norm_∂x_c², c)
    + ▶y_aca(i, j, k, grid, norm_∂y_c², c)
    + ▶z_aac(i, j, k, grid, norm_∂z_c², c)
)

"""
    ∇_κ_∇T(i, j, k, grid, T, closure, diffusivities)

Return the diffusive flux divergence `∇ ⋅ (κ ∇ c)` for the turbulence
`closure`, where `c` is an array of scalar data located at cell centers.
"""
@inline ∇_κ_∇T(i, j, k, grid, T, closure::VAMD, diffusivities) = (
      ∂x_caa(i, j, k, grid, κ_∂x_c, T, diffusivities.κₑ.T, closure)
    + ∂y_aca(i, j, k, grid, κ_∂y_c, T, diffusivities.κₑ.T, closure)
    + ∂z_aac(i, j, k, grid, κ_∂z_c, T, diffusivities.κₑ.T, closure)
)

"""
    ∇_κ_∇S(i, j, k, grid, S, closure, diffusivities)

Return the diffusive flux divergence `∇ ⋅ (κ ∇ S)` for the turbulence
`closure`, where `c` is an array of scalar data located at cell centers.
"""
@inline ∇_κ_∇S(i, j, k, grid, S, closure::VAMD, diffusivities) = (
      ∂x_caa(i, j, k, grid, κ_∂x_c, S, diffusivities.κₑ.S, closure)
    + ∂y_aca(i, j, k, grid, κ_∂y_c, S, diffusivities.κₑ.S, closure)
    + ∂z_aac(i, j, k, grid, κ_∂z_c, S, diffusivities.κₑ.S, closure)
)

function calc_diffusivities!(K, grid, closure::VAMD, buoyancy, U, Φ)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds K.νₑ[i, j, k]   = ν_ccc(i, j, k, grid, closure, nothing, buoyancy, U, Φ)
                @inbounds K.κₑ.T[i, j, k] = κ_ccc(i, j, k, grid, closure, Φ.T,     buoyancy, U, Φ)
                @inbounds K.κₑ.S[i, j, k] = κ_ccc(i, j, k, grid, closure, Φ.S,     buoyancy, U, Φ)
            end
        end
    end
    return nothing
end

#####
##### Filter width at various locations
#####

# Recall that filter widths are 2x the grid spacing in VAMD
@inline Δᶠx_ccc(i, j, k, grid::RegularCartesianGrid) = 2 * grid.Δx
@inline Δᶠy_ccc(i, j, k, grid::RegularCartesianGrid) = 2 * grid.Δy
@inline Δᶠz_ccc(i, j, k, grid::RegularCartesianGrid) = 2 * grid.Δz

for loc in (:ccf, :fcc, :cfc, :ffc, :cff, :fcf), ξ in (:x, :y, :z)
    Δ_loc = Symbol(:Δᶠ, ξ, :_, loc)
    Δ_ccc = Symbol(:Δᶠ, ξ, :_ccc)
    @eval begin
        const $Δ_loc = $Δ_ccc
    end
end


text:
load("cuplus.jl")
load("cublas.jl")
load("cufft.jl")

max_stars_repo_name: stefan-k/julia-CuMatrix

max_issues_repo_name: stefan-k/julia-CuMatrix

max_forks_repo_name: stefan-k/julia-CuMatrix