gistfile1.txt

function build_full_finite_difference_matrix()
    result = spzeros(Int64,64,64)
    function to1d(i,j,k)
        return i + j*4 + k*4*4
    end

    #Point values
    for i in 0:1, j in 0:1, k in 0:1
        aindex::Int64 = i + 2*j + 4*k +1
        result[aindex,to1d(i+1,j+1,k+1)+1] = 8
    end


    #Specify first derivatives
    for i in 0:1, j in 0:1, k in 0:1
        for ddirection in 1:3
            if ddirection == 1
                ddirectionT::Tuple{Int64,Int64,Int64} = (1,0,0)
            elseif ddirection == 2
                ddirectionT = (0,1,0)
            elseif ddirection == 3
                ddirectionT = (0,0,1)
            end
            #res += Us[earthIndex(current_index .+ ddirectionT,nx,ny,nt)]
            #res -= Us[earthIndex(current_index .- ddirectionT,nx,ny,nt)]
            aindex::Int64 = i + 2*j + 4*k + 8*ddirection + 1
            result[aindex, to1d(
                i + ddirectionT[1] + 1, j + ddirectionT[2] + 1, k + ddirectionT[3] + 1)+1
                ] = 4
            result[aindex,to1d(
                i - ddirectionT[1] + 1, j - ddirectionT[2] + 1, k + ddirectionT[3] + 1)+1
            ] = -4
        end
    end

    for i in 0:1, j in 0:1, k in 0:1
        for ddirection1 in 1:2
            if ddirection1 == 1
                ddirection1T::Tuple{Int64,Int64,Int64} = (1,0,0)
            else
                ddirection1T = (0,1,0)
            end
            for ddirection2 in (ddirection1+1):3
                if ddirection2 == 2
                    ddirection2T::Tuple{Int64,Int64,Int64} = (0,1,0)
                else
                    ddirection2T = (0,0,1)
                end

                    aindex = i + 2*j + 4*k + (2*(ddirection1-1) + (ddirection2 - ddirection1 - 1) + 4)*8 + 1
                    result[aindex,to1d(i + ddirection1T[1] + ddirection2T[1] + 1,
                             j + ddirection1T[2] + ddirection2T[2] + 1,
                             k + ddirection1T[3] + ddirection2T[3] + 1
                             ) + 1 ] = 2

                    result[aindex,to1d(i + ddirection1T[1] - ddirection2T[1] + 1,
                            j + ddirection1T[2] - ddirection2T[2] + 1,
                            k + ddirection1T[3] - ddirection2T[3] + 1
                            ) + 1] = 2

                    result[aindex,to1d(i - ddirection1T[1] + ddirection2T[1] + 1,
                            j - ddirection1T[2] + ddirection2T[2] + 1,
                            k - ddirection1T[3] + ddirection2T[3] + 1
                            ) + 1] = -2

                    result[aindex,to1d(i - ddirection1T[1] - ddirection2T[1] + 1,
                            j - ddirection1T[2] - ddirection2T[2] + 1,
                            k - ddirection1T[3] - ddirection2T[3] + 1
                            ) + 1] = 2
            end
        end
    end
    #Specfiy (mixed) third derivatives
    for i in 0:1, j in 0:1, k in 0:1
        aindex = i + 2*j + 4*k + 56 + 1

        res = 0.0
        result[aindex,to1d(i+1+1,j+1+1,k+1+1) + 1] = 1
        result[aindex,to1d(i+1+1,j+1+1,k-1+1) + 1] = -1
        result[aindex,to1d(i+1+1,j-1+1,k+1+1) + 1] = -1
        result[aindex,to1d(i+1+1,j-1+1,k-1+1) + 1] = 1
        result[aindex,to1d(i-1+1,j+1+1,k+1+1) + 1] = -1
        result[aindex,to1d(i-1+1,j+1+1,k-1+1) + 1] = 1
        result[aindex,to1d(i-1+1,j-1+1,k+1+1) + 1] = 1
        result[aindex,to1d(i-1+1,j-1+1,k-1+1) + 1] = -1
    end
    return result
end