module test1
import LinearAlgebra.det
import LinearAlgebra.I
eye(n) = Matrix{Float64}(I,n,n)


struct BulkMesh
    nodepos::Array{Float64, 2}
    eltnodes::Array{Int, 2}
end

struct BulkMeshAux
    B_matrix::Vector{Vector{Array{Float64,2}}}
end

function problemsetup0(displacement_delta)
    spacedim = 2
    meshfilename = "testmesh.dat"
    nnode0, nelt0, nodes0, elements0, labels = readmesh(meshfilename)
    nodes0_sample = Vector{Vector{Float64}}()
    nodes0_pin1 = Vector{Vector{Float64}}()
    nodes0_pin3 = Vector{Vector{Float64}}()
    node0lab = zeros(Int, nnode0)
    node0renum = zeros(Int, nnode0)

    for j = 1 : nnode0
        push!(nodes0_sample, nodes0[j])
        node0lab[j] = 10
        node0renum[j] = length(nodes0_sample)

    end

    elements0_sample = Vector{Vector{Int}}()
    elements0_pin1 = Vector{Vector{Int}}()
    elements0_pin3 = Vector{Vector{Int}}()

    for j = 1 : nelt0
        erenum = node0renum[elements0[j]]
        push!(elements0_sample, erenum)
    end

    bulkmesh0_sample = make_bulkmesh0(nodes0_sample,
                                      elements0_sample)
    nnode0_sample = size(bulkmesh0_sample.nodepos, 1)
    bulkmesh0_pin1 = make_bulkmesh0(nodes0_pin1,
                                    elements0_pin1)
    nnode_pin1 = size(bulkmesh0_pin1.nodepos, 1)
    bulkmesh0_pin3 = make_bulkmesh0(nodes0_pin3,
                                    elements0_pin3)
    bulkmesh_sample = deepcopy(bulkmesh0_sample)
    nnode_sample = size(bulkmesh_sample.nodepos,1)
    allnodes = vcat(bulkmesh_sample.nodepos,
                    bulkmesh0_pin1.nodepos,
                    bulkmesh0_pin3.nodepos)
    allelts = vcat(bulkmesh_sample.eltnodes,
                   bulkmesh0_pin1.eltnodes .+ nnode_sample,
                   bulkmesh0_pin3.eltnodes .+ (nnode_sample + nnode_pin1))

    completemesh = BulkMesh(allnodes,
                            allelts)
    bulkmesh_aux = make_bulkmesh_aux(completemesh, triquad3pt())
    numstep = Int(ceil(2.5 / displacement_delta))
    if nnode0_sample == 3
        println("nnode0_sample has changed")
        error("W")
    end

    completemesh
end

function readmesh(filename)
    nnode0 = 0
    nelt0 = 0
    nodes = Vector{Vector{Float64}}()
    elements = Vector{Vector{Int}}()
    labels = Vector{Int}()
    open(filename, "r") do readhandle
        s = readline(readhandle)
        w = split(s)
        spacedim = parse(Int, w[1])
        nnode_per_elt = parse(Int, w[2])
        s = readline(readhandle)
        w = split(s)
        nnode0 = parse(Int, w[1])
        nelt0 = parse(Int, w[2])
        for nn = 1 : nnode0
            s = readline(readhandle)
            w = split(s)
            push!(nodes, zeros(0))
            for d = 1 : spacedim
                push!(nodes[nn], parse(Float64, w[d]))
            end
            push!(labels, parse(Int, w[spacedim+1]))
        end
        for ne = 1 : nelt0
            s = readline(readhandle)
            w = split(s)
            push!(elements, zeros(Int,0))
            for i = 1 : nnode_per_elt
                push!(elements[ne], parse(Int, w[i]))
            end
        end
    end
    return nnode0, nelt0, nodes, elements, labels
end

function make_bulkmesh_aux(bulkmesh0, quadrule)
    spacedim = 2
    xyr, polyv1, shapefunc, shapefuncdxi1, shapefuncdxi2 = t6_shapefunc()
    ngau2 = size(quadrule.points,1)
    NIxidxi = Vector{Array{Float64,2}}()
    for s = 1 : ngau2
        qpxi = quadrule.points[s,1]
        qpeta = quadrule.points[s,2]
        monomialvec = [1.0 qpxi  qpeta]
        push!(NIxidxi,[monomialvec * shapefuncdxi1
                       monomialvec * shapefuncdxi2])
    end
    nelement = size(bulkmesh0.eltnodes,1)
    B_matrix = Vector{Vector{Array{Float64,2}}}()
    eye2 = eye(spacedim)
    for e = 1 : nelement
        push!(B_matrix, Vector{Array{Float64,2}}())
        mynodes = bulkmesh0.nodepos[vec(bulkmesh0.eltnodes[e,:]),:]
        for s = 1 : ngau2
            qpxi = quadrule.points[s,1]
            qpeta = quadrule.points[s,2]
            monomialvec = [1 qpxi qpeta]
            dXdxi = [monomialvec * shapefuncdxi1 * mynodes
                     monomialvec * shapefuncdxi2 * mynodes]'
            invdXdxi = inv(dXdxi)
            detdXdxi = det(dXdxi)
            B = invdXdxi' * NIxidxi[s]
            push!(B_matrix[e], B)
        end
    end
    bulkmesh_aux = BulkMeshAux(B_matrix)
    return bulkmesh_aux
end


function t6_shapefunc()
    xyr = [0.0 0.0
           0.5 0.0
           1.0 0.0
           0.5 0.5
           0.0 1.0
           0.0 0.5]
    polyv1 = [ones(6,1) xyr[:,1] xyr[:,2] xyr[:,1].^2  xyr[:,1].*xyr[:,2] xyr[:,2].^2]
    shapefunc = inv(polyv1)
    shapefuncdxi1 = vcat(vec(shapefunc[2,:])', 2.0*vec(shapefunc[4,:])', vec(shapefunc[5,:])')
    shapefuncdxi2 = vcat(vec(shapefunc[3,:])', vec(shapefunc[5,:])', 2.0*vec(shapefunc[6,:])')
    xyr, polyv1, shapefunc, shapefuncdxi1, shapefuncdxi2
end


function make_bulkmesh0(nodes0,
                        elements0)
    nnode0 = length(nodes0)
    nelt0 = length(elements0)
    spacedim  = 2
    nnode_per_element = 6
    bulkmesh0 = BulkMesh(zeros(nnode0, spacedim),
                         zeros(Int, nelt0, nnode_per_element))

    for nn = 1 : nnode0
        bulkmesh0.nodepos[nn,:] .= nodes0[nn][:]
    end
    for ne = 1 : nelt0
        bulkmesh0.eltnodes[ne,:] .= elements0[ne][:]
    end
    bulkmesh0
end

struct QuadRule
    points::Array{Float64,2}
    weights::Vector{Float64}
end


# 3-point quadrature for T6 elements

function triquad3pt()
    QuadRule([1.0/6.0 1.0/6.0; 1.0/6.0 2.0/3.0; 2.0/3.0 1.0/6.0],
             [1.0/6.0,  1.0/6.0,  1.0/6.0])
end


end

for k = 1 : 500
    bulkmesh = test1.problemsetup0(.1)
    println(bulkmesh.nodepos[1,1])
end