text:
module hughes_cantilever_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.AlgoBaseModule: evalconvergencestudy
using FinEtoolsDeforLinear.AlgoDeforLinearModule: linearstatics, exportstresselementwise, exportstress
using Statistics: mean
using LinearAlgebra: Symmetric, cholesky

# Example from TJR Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, 1987.
# Two-dimensional plane-elasticity Solution of beam deflection. Shear force P at one end, and corresponding reactions (tractions) at the other end. The Analytical solution is plane-strain.
# Find out: where does the exact solution come from?

# Isotropic material
E=1.0;
P=1.0;
L=16.0;
c=2.0;
h=1.0;
W=2/3*h*c^3;
CTE = 0.0
tolerance = 0.00001*c

function getfrcL!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    copyto!(forceout, [0.0; 0.0; P/2/W*(c^2-XYZ[3]^2)])
end

function getfrc0!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    copyto!(forceout, [P*L/W*XYZ[3]; 0.0; -P/2/W*(c^2-XYZ[3]^2)])
end

function evaluateerrors(filebase, modeldatasequence)
    println("")
    println("Stress RMS error")
    for md = modeldatasequence
        md["targetfields"] = [e["field"] for e in md["postprocessing"]["exported"]]
    end
    elementsizes, errornorms, p = evalconvergencestudy(modeldatasequence)

    println("Normalized Approximate Error = $(errornorms)")

    f = log.(vec(errornorms))
    A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
    p = A \ f
    println("Linear log-log fit: p = $(p)")

    csvFile = filebase * "_Stress" * ".CSV"
    savecsv(csvFile,
    elementsizes=vec(elementsizes[1:end-1]),
    elementsizes2=vec(elementsizes[1:end-1].^2),
    elementsizes3=vec(elementsizes[1:end-1].^3),
    errornorms=vec(errornorms)
    )
    println("Wrote $csvFile")

    println("")
    println("Displacement RMS error")
    for md = modeldatasequence
        md["targetfields"] = [md["u"] for r in md["regions"]]
    end
    elementsizes, errornorms, p = evalconvergencestudy(modeldatasequence)

    println("Normalized Approximate Error = $(errornorms)")

    f = log.(vec(errornorms))
    A = hcat(log.(vec(elementsizes[1:end-1])), ones(size(f)))
    p = A \ f
    println("Linear log-log fit: p = $(p)")

    csvFile = filebase * "_Displ" * ".CSV"
    savecsv(csvFile,
    elementsizes=vec(elementsizes[1:end-1]),
    elementsizes2=vec(elementsizes[1:end-1].^2),
    elementsizes3=vec(elementsizes[1:end-1].^3),
    errornorms=vec(errornorms)
    )
    println("Wrote $csvFile")
end

function hughes_cantilever_stresses_H8_by_hand()
    elementtag = "H8"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)
    nu=0.3; # COMPRESSIBLE
    E1=E/(1-nu^2);  # Plane strain
    nu1=nu/(1-nu); I=h*(2*c)^3/12;
    exactux(x,y) = (P/(6*E1*I)*(-y)*(3*(L^2-(L-x)^2)+(2+nu1)*(y^2-c^2)));
    exactuy(x,y) = (P/(6*E1*I)*(((L-x)^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x)*y^2));

    n = 2 #

    nL = 3*n # number of elements lengthwise
    nc = 2*n # number of elements through the depth
    nh = 1 # number of elements through the thickness
    xs = collect(linearspace(0.0, L, nL+1))
    ys = collect(linearspace(0.0, h, nh+1))
    zs = collect(linearspace(-c, +c, nc+1))
    fens,fes = H8blockx(xs, ys, zs)
    # fens,fes = H8toH20(fens,fes)
    bfes = meshboundary(fes)
    # end cross-section surface  for the shear loading
    sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
    # 0 cross-section surface  for the reactions
    sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, 0.0, E, nu, CTE)

    # Material orientation matrix
    csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(csmatout, csmat)
    end

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    setebc!(u,lx0,true,1,0.0)
    setebc!(u,lx0,true,2,0.0)
    setebc!(u,lx0,true,3,0.0)
    lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
    lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
    setebc!(u,vcat(lx1, lx2),true,1,0.0)
    ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
    setebc!(u,vcat(ly1, ly2),true,2,0.0)
    applyebc!(u)
    numberdofs!(u)

    fi = ForceIntensity(FFlt, 3, getfrc0!);
    el1femm = FEMMBase(IntegDomain(subset(bfes, sshear0), GaussRule(2, 2)))
    F1 = distribloads(el1femm, geom, u, fi, 2);
    fi = ForceIntensity(FFlt, 3, getfrcL!);
    el2femm = FEMMBase(IntegDomain(subset(bfes, sshearL), GaussRule(2, 2)))
    F2 = distribloads(el2femm, geom, u, fi, 2);

    associategeometry!(femm, geom)
    K = stiffness(femm, geom, u)

    K = cholesky(K)
    U = K\(F1 + F2)
    scattersysvec!(u,U[:])

    Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
    utip = mean(u.values[Tipl, 3])
    println("Deflection: $(utip), compared to $(exactuy(L,0.0))")

    File =  "hughes_cantilever_stresses_H8_by_hand.vtk"
    vtkexportmesh(File, fens, fes;  vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)
    # modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)", "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>[5])
    # modeldata = exportstresselementwise(modeldata)

    # modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
    # "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
    # "component"=>collect(1:6))
    # modeldata = exportstresselementwise(modeldata)
    # stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

    true

end # hughes_cantilever_stresses_H8_by_hand

function hughes_cantilever_stresses_H20_by_hand()
    elementtag = "H20"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)
    nu=0.3; # COMPRESSIBLE
    E1=E/(1-nu^2);  # Plane strain
    nu1=nu/(1-nu); I=h*(2*c)^3/12;
    exactux(x,y) = (P/(6*E1*I)*(-y)*(3*(L^2-(L-x)^2)+(2+nu1)*(y^2-c^2)));
    exactuy(x,y) = (P/(6*E1*I)*(((L-x)^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x)*y^2));

    n = 8 #

    nL = 3*n # number of elements lengthwise
    nc = 2*n # number of elements through the depth
    nh = 1 # number of elements through the thickness
    xs = collect(linearspace(0.0, L, nL+1))
    ys = collect(linearspace(0.0, h, nh+1))
    zs = collect(linearspace(-c, +c, nc+1))
    fens,fes = H8blockx(xs, ys, zs)
    fens,fes = H8toH20(fens,fes)
    bfes = meshboundary(fes)
    # end cross-section surface  for the shear loading
    sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
    # 0 cross-section surface  for the reactions
    sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, 0.0, E, nu, CTE)

    # Material orientation matrix
    csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(csmatout, csmat)
    end

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    setebc!(u,lx0,true,1,0.0)
    setebc!(u,lx0,true,2,0.0)
    setebc!(u,lx0,true,3,0.0)
    lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
    lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
    setebc!(u,vcat(lx1, lx2),true,1,0.0)
    ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
    setebc!(u,vcat(ly1, ly2),true,2,0.0)
    applyebc!(u)
    numberdofs!(u)
    u

    fi = ForceIntensity(FFlt, 3, getfrc0!);
    el1femm = FEMMBase(IntegDomain(subset(bfes, sshear0), GaussRule(2, 2)))
    F1 = distribloads(el1femm, geom, u, fi, 2);
    fi = ForceIntensity(FFlt, 3, getfrcL!);
    el2femm = FEMMBase(IntegDomain(subset(bfes, sshearL), GaussRule(2, 2)))
    F2 = distribloads(el2femm, geom, u, fi, 2);

    associategeometry!(femm, geom)
    K = stiffness(femm, geom, u)

    K = cholesky(K)
    U = K\(F1 + F2)
    scattersysvec!(u,U[:])

    Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
    utip = mean(u.values[Tipl, 3])
    println("Deflection: $(utip), compared to $(exactuy(L,0.0))")

    File =  "hughes_cantilever_stresses_H8_by_hand.vtk"
    vtkexportmesh(File, fens, fes;  vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)
    # modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)", "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>[5])
    # modeldata = exportstresselementwise(modeldata)

    # modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
    # "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
    # "component"=>collect(1:6))
    # modeldata = exportstresselementwise(modeldata)
    # stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

    true

end # hughes_cantilever_stresses_MST10

function hughes_cantilever_stresses_T10_by_hand()
    elementtag = "T10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)
    nu=0.3; # COMPRESSIBLE
    E1=E/(1-nu^2);  # Plane strain
    nu1=nu/(1-nu); I=h*(2*c)^3/12;
    exactux(x,y) = (P/(6*E1*I)*(-y)*(3*(L^2-(L-x)^2)+(2+nu1)*(y^2-c^2)));
    exactuy(x,y) = (P/(6*E1*I)*(((L-x)^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x)*y^2));

    n = 2 #

    nL = 3*n # number of elements lengthwise
    nc = 2*n # number of elements through the depth
    nh = n # number of elements through the thickness
    xs = collect(linearspace(0.0, L, nL+1))
    ys = collect(linearspace(0.0, h, nh+1))
    zs = collect(linearspace(-c, +c, nc+1))
    fens,fes = T10blockx(xs, ys, zs)
    bfes = meshboundary(fes)
    # end cross-section surface  for the shear loading
    sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
    # 0 cross-section surface  for the reactions
    sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, 0.0, E, nu, CTE)

    # Material orientation matrix
    csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(csmatout, csmat)
    end

    gr = SimplexRule(3, 4)

    femm = FEMMDeforLinear(MR, IntegDomain(fes, gr), material)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    setebc!(u,lx0,true,1,0.0)
    setebc!(u,lx0,true,2,0.0)
    setebc!(u,lx0,true,3,0.0)
    lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
    lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
    setebc!(u,vcat(lx1, lx2),true,1,0.0)
    ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
    setebc!(u,vcat(ly1, ly2),true,2,0.0)
    applyebc!(u)
    numberdofs!(u)

    fi = ForceIntensity(FFlt, 3, getfrc0!);
    el1femm = FEMMBase(IntegDomain(subset(bfes, sshear0), SimplexRule(2, 3)))
    F1 = distribloads(el1femm, geom, u, fi, 2);
    fi = ForceIntensity(FFlt, 3, getfrcL!);
    el2femm = FEMMBase(IntegDomain(subset(bfes, sshearL), SimplexRule(2, 3)))
    F2 = distribloads(el2femm, geom, u, fi, 2);

    associategeometry!(femm, geom)
    K = stiffness(femm, geom, u)
    K=cholesky(K)
    U = K\(F1 + F2)
    scattersysvec!(u,U[:])

    Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
    utip = mean(u.values[Tipl, 3])
    println("Deflection: $(utip), compared to $(exactuy(L,0.0))")

    File =  "hughes_cantilever_stresses_T10_by_hand.vtk"
    vtkexportmesh(File, fens, fes;  vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)
    # modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)", "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>[5])
    # modeldata = exportstresselementwise(modeldata)

    # modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
    # "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
    # "component"=>collect(1:6))
    # modeldata = exportstresselementwise(modeldata)
    # stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

    true

end # hughes_cantilever_stresses_MST10

function hughes_cantilever_stresses_MST10_by_hand()
    elementtag = "MST10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)
    nu=0.3; # COMPRESSIBLE
    E1=E/(1-nu^2);  # Plane strain
    nu1=nu/(1-nu); I=h*(2*c)^3/12;
    exactux(x,y) = (P/(6*E1*I)*(-y)*(3*(L^2-(L-x)^2)+(2+nu1)*(y^2-c^2)));
    exactuy(x,y) = (P/(6*E1*I)*(((L-x)^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x)*y^2));

    n = 2 #

    nL = 3*n # number of elements lengthwise
    nc = 2*n # number of elements through the depth
    nh = n # number of elements through the thickness
    xs = collect(linearspace(0.0, L, nL+1))
    ys = collect(linearspace(0.0, h, nh+1))
    zs = collect(linearspace(-c, +c, nc+1))
    fens,fes = T10blockx(xs, ys, zs)
    bfes = meshboundary(fes)
    # end cross-section surface  for the shear loading
    sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
    # 0 cross-section surface  for the reactions
    sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, 0.0, E, nu, CTE)

    # Material orientation matrix
    csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copyto!(csmatout, csmat)
    end

    gr = SimplexRule(3, 4)

    femm = FEMMDeforLinearMST10(MR, IntegDomain(fes, gr), material)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    setebc!(u,lx0,true,1,0.0)
    setebc!(u,lx0,true,2,0.0)
    setebc!(u,lx0,true,3,0.0)
    lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
    lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
    setebc!(u,vcat(lx1, lx2),true,1,0.0)
    ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
    setebc!(u,vcat(ly1, ly2),true,2,0.0)
    applyebc!(u)
    numberdofs!(u)

    fi = ForceIntensity(FFlt, 3, getfrc0!);
    el1femm = FEMMBase(IntegDomain(subset(bfes, sshear0), SimplexRule(2, 3)))
    F1 = distribloads(el1femm, geom, u, fi, 2);
    fi = ForceIntensity(FFlt, 3, getfrcL!);
    el2femm = FEMMBase(IntegDomain(subset(bfes, sshearL), SimplexRule(2, 3)))
    F2 = distribloads(el2femm, geom, u, fi, 2);

    associategeometry!(femm, geom)
    K =stiffness(femm, geom, u)
    K=cholesky(K)
    U=  K\(F1 + F2)
    scattersysvec!(u,U[:])

    Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
    utip = mean(u.values[Tipl, 3])
    println("Deflection: $(utip), compared to $(exactuy(L,0.0))")

    # modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)", "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>[5])
    # modeldata = exportstresselementwise(modeldata)

    # modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
    # "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
    # "component"=>collect(1:6))
    # modeldata = exportstresselementwise(modeldata)
    # stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

    true

end # hughes_cantilever_stresses_MST10

function hughes_cantilever_stresses_MST10()
    elementtag = "MST10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)
    nu=0.3; # COMPRESSIBLE
    E1=E/(1-nu^2);  # Plane strain
    nu1=nu/(1-nu); I=h*(2*c)^3/12;
    exactux(x,y) = (P/(6*E1*I)*(-y)*(3*(L^2-(L-x)^2)+(2+nu1)*(y^2-c^2)));
    exactuy(x,y) = (P/(6*E1*I)*(((L-x)^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x)*y^2));

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the depth
        nh = n # number of elements through the thickness
        xs = collect(linearspace(0.0, L, nL+1))
        ys = collect(linearspace(0.0, h, nh+1))
        zs = collect(linearspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR, 0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinearMST10(MR, IntegDomain(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        Trac0 = FDataDict("traction_vector"=>getfrc0!, "femm"=>FEMMBase(IntegDomain(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!, "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip), compared to $(exactuy(L,0.0))")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)", "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy, "component"=>[5])
        modeldata = exportstresselementwise(modeldata)

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = exportstresselementwise(modeldata)
        stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_MST10

function hughes_cantilever_stresses_MST10_incompressible()
    elementtag = "MST10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)
    nu=0.499999999; # INCOMPRESSIBLE
    E1=E/(1-nu^2);  # Plane strain
    nu1=nu/(1-nu); I=h*(2*c)^3/12;
    exactux(x,y) = (P/(6*E1*I)*(-y)*(3*(L^2-(L-x)^2)+(2+nu1)*(y^2-c^2)));
    exactuy(x,y) = (P/(6*E1*I)*(((L-x)^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x)*y^2));

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the wwith
        nh = n # number of elements through the thickness
        xs = collect(linearspace(0.0, L, nL+1))
        ys = collect(linearspace(0.0, h, nh+1))
        zs = collect(linearspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR,  0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinearMST10(MR,
        IntegDomain(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        Trac0 = FDataDict("traction_vector"=>getfrc0!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip), compared to $(exactuy(L,0.0))")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>[5])
        modeldata = exportstresselementwise(modeldata)

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = exportstresselementwise(modeldata)
        stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_incompressible_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_MST10_incompressible

function hughes_cantilever_stresses_nodal_MST10()
    elementtag = "MST10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)
    nu=0.3; # COMPRESSIBLE
    E1=E/(1-nu^2);  # Plane strain
    nu1=nu/(1-nu); I=h*(2*c)^3/12;
    exactux(x,y) = (P/(6*E1*I)*(-y)*(3*(L^2-(L-x)^2)+(2+nu1)*(y^2-c^2)));
    exactuy(x,y) = (P/(6*E1*I)*(((L-x)^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x)*y^2));

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the wwith
        nh = n # number of elements through the thickness
        xs = collect(linearspace(0.0, L, nL+1))
        ys = collect(linearspace(0.0, h, nh+1))
        zs = collect(linearspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR,   0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4) # rule for tetrahedral meshes

        region = FDataDict("femm"=>FEMMDeforLinearMST10(MR,
        IntegDomain(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        Trac0 = FDataDict("traction_vector"=>getfrc0!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip), compared to $(exactuy(L,0.0))")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_nodal_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6), "nodevalmethod"=>:averaging, "reportat"=>:extraptrend)
        modeldata = exportstress(modeldata)
        stressfields = NodalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_nodal_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_nodal_MST10

function hughes_cantilever_stresses_nodal_T10()
    elementtag = "T10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)
    nu=0.3; # COMPRESSIBLE
    E1=E/(1-nu^2);  # Plane strain
    nu1=nu/(1-nu); I=h*(2*c)^3/12;
    exactux(x,y) = (P/(6*E1*I)*(-y)*(3*(L^2-(L-x)^2)+(2+nu1)*(y^2-c^2)));
    exactuy(x,y) = (P/(6*E1*I)*(((L-x)^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x)*y^2));

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the wwith
        nh = n # number of elements through the thickness
        xs = collect(linearspace(0.0, L, nL+1))
        ys = collect(linearspace(0.0, h, nh+1))
        zs = collect(linearspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR,
        0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? 1.0 : 0.0 for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinear(MR,
        IntegDomain(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        Trac0 = FDataDict("traction_vector"=>getfrc0!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip), compared to $(exactuy(L,0.0))")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_nodal_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = exportstress(modeldata)
        stressfields = NodalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_nodal_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_nodal_T10

function hughes_cantilever_stresses_T10()
    elementtag = "T10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)
    nu=0.3; # COMPRESSIBLE
    E1=E/(1-nu^2);  # Plane strain
    nu1=nu/(1-nu); I=h*(2*c)^3/12;
    exactux(x,y) = (P/(6*E1*I)*(-y)*(3*(L^2-(L-x)^2)+(2+nu1)*(y^2-c^2)));
    exactuy(x,y) = (P/(6*E1*I)*(((L-x)^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x)*y^2));

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the wwith
        nh = n # number of elements through the thickness
        xs = collect(linearspace(0.0, L, nL+1))
        ys = collect(linearspace(0.0, h, nh+1))
        zs = collect(linearspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR,
        0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        Trac0 = FDataDict("traction_vector"=>getfrc0!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip), compared to $(exactuy(L,0.0))")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>[5])
        modeldata = exportstresselementwise(modeldata)

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = exportstresselementwise(modeldata)
        stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_T10

function hughes_cantilever_stresses_T10_incompressible()
    elementtag = "T10"
    println("""
    Cantilever example.  Hughes 1987. Element: $(elementtag)
    """)
        nu=0.499999999; # INCOMPRESSIBLE
        E1=E/(1-nu^2);  # Plane strain
        nu1=nu/(1-nu); I=h*(2*c)^3/12;
        exactux(x,y) = (P/(6*E1*I)*(-y)*(3*(L^2-(L-x)^2)+(2+nu1)*(y^2-c^2)));
        exactuy(x,y) = (P/(6*E1*I)*(((L-x)^3-L^3)-((4+5*nu1)*c^2+3*L^2)*(L-x-L)+3*nu1*(L-x)*y^2));

    modeldatasequence = FDataDict[]
    for n = [1 2 4 8] #

        nL = 3*n # number of elements lengthwise
        nc = 2*n # number of elements through the wwith
        nh = n # number of elements through the thickness
        xs = collect(linearspace(0.0, L, nL+1))
        ys = collect(linearspace(0.0, h, nh+1))
        zs = collect(linearspace(-c, +c, nc+1))
        fens,fes = T10blockx(xs, ys, zs)
        bfes = meshboundary(fes)
        # end cross-section surface  for the shear loading
        sshearL = selectelem(fens, bfes; facing=true, direction = [+1.0 0.0 0.0])
        # 0 cross-section surface  for the reactions
        sshear0 = selectelem(fens, bfes; facing=true, direction = [-1.0 0.0 0.0])

        MR = DeforModelRed3D
        material = MatDeforElastIso(MR,  0.0, E, nu, CTE)

        # Material orientation matrix
        csmat = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            copyto!(csmatout, csmat)
        end

        gr = SimplexRule(3, 4)

        region = FDataDict("femm"=>FEMMDeforLinear(MR,
        IntegDomain(fes, gr), material))

        lx0 = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
        # println("lx0 = $(lx0)")
        ex01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
        ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
        ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
        lx1 = selectnode(fens, box=[0.0 0.0 0.0 0.0 c c], inflate=tolerance)
        lx2 = selectnode(fens, box=[0.0 0.0 0.0 0.0 -c -c], inflate=tolerance)
        # println("vcat(lx1, lx2) = $(vcat(lx1, lx2))")
        ex04 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>vcat(lx1, lx2) )
        ly1 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        ly2 = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        # println("vcat(ly1, ly2) = $(vcat(ly1, ly2))")
        ey01 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly1, ly2) )

        Trac0 = FDataDict("traction_vector"=>getfrc0!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshear0), SimplexRule(2, 3))))
        TracL = FDataDict("traction_vector"=>getfrcL!,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, sshearL), SimplexRule(2, 3))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>[region],
        "essential_bcs"=>[ex01, ex02, ex03, ex04, ey01],
        "traction_bcs"=>[Trac0, TracL],
        "temperature_change"=>FDataDict("temperature"=>0.0)
        )
        modeldata = linearstatics(modeldata)

        u = modeldata["u"]
        geom = modeldata["geom"]

        Tipl = selectnode(fens, box=[L L 0.0 0.0 0. 0.], inflate=tolerance)
        utip = mean(u.values[Tipl, 3])
        println("Deflection: $(utip), compared to $(exactuy(L,0.0))")

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>[5])
        modeldata = exportstresselementwise(modeldata)

        modeldata["postprocessing"] = FDataDict("file"=>"hughes_cantilever_stresses_$(elementtag)",
        "outputcsys"=>CSys(3, 3, updatecs!), "quantity"=>:Cauchy,
        "component"=>collect(1:6))
        modeldata = exportstresselementwise(modeldata)
        stressfields = ElementalField[modeldata["postprocessing"]["exported"][1]["field"]]

        modeldata["geometricaltolerance"] = tolerance
        modeldata["elementsize"] = 1.0 / n
        push!(modeldatasequence, modeldata)
    end

    filebase = "hughes_cantilever_stresses_incompressible_$(elementtag)"
    evaluateerrors(filebase, modeldatasequence)

    true

end # hughes_cantilever_stresses_T10_incompressible

function allrun()
    println("#####################################################")
    println("# hughes_cantilever_stresses_H8_by_hand ")
    hughes_cantilever_stresses_H8_by_hand()
    println("#####################################################")
    println("# hughes_cantilever_stresses_H20_by_hand ")
    hughes_cantilever_stresses_H20_by_hand()
    println("#####################################################")
    println("# hughes_cantilever_stresses_T10_by_hand ")
    hughes_cantilever_stresses_T10_by_hand()
    println("#####################################################")
    println("# hughes_cantilever_stresses_MST10 ")
    hughes_cantilever_stresses_MST10()
    println("#####################################################")
    println("# hughes_cantilever_stresses_MST10_incompressible ")
    hughes_cantilever_stresses_MST10_incompressible()
    println("#####################################################")
    println("# hughes_cantilever_stresses_nodal_MST10 ")
    hughes_cantilever_stresses_nodal_MST10()
    println("#####################################################")
    println("# hughes_cantilever_stresses_nodal_T10 ")
    hughes_cantilever_stresses_nodal_T10()
    println("#####################################################")
    println("# hughes_cantilever_stresses_T10 ")
    hughes_cantilever_stresses_T10()
    println("#####################################################")
    println("# hughes_cantilever_stresses_T10_incompressible ")
    hughes_cantilever_stresses_T10_incompressible()
    return true
end # function allrun

end # module hughes_cantilever_examples


text:
module Parser

using CombinedParsers
using TextParse

# @syntax subterm = Either{Rational{Int}}(Any[TextParse.Numeric(Int)]);
# @syntax for parenthesis in subterm
#     mult         = evaluate |> join(subterm, CharIn("*/"), infix=:prefix )
#     @syntax term = evaluate |> join(mult,    CharIn("+-"), infix=:prefix )
#     Sequence(2,'(',term,')')
# end;

end

ext: jl

max_stars_repo_name: mesa-x/spreadsheet

max_issues_repo_name: mesa-x/spreadsheet

max_forks_repo_name: mesa-x/spreadsheet


text:
# Example 1 of Tube Based EMPC Paper
#TODO clean up and add compute_tightened_sets to library!

using MPCHelper

using JuMP
using Plots
using Ipopt
using LinearAlgebra
using LazySets

## Initialize System
A = [0.7776 -0.0045; 26.6185 1.8555]
B = [-0.0004 0.2907]'
W = [-0.0002 0.0893; 0.139 1.2267]
K = [-0.6457 -5.4157]

Aʷ = [1 0.0; 0 1; -1 0; 0 -1]
bʷ = [0.1; 0.1; 0.1; 0.1]
 = linear_map(W, HPolyhedron(Aʷ, bʷ)) # corresponds to W⁻¹Aʷ <= bʷ
plot(𝒲, 1e-3, alpha=0.1)
plot!(𝒲, 1e-3, alpha=0.1)
a = 𝒲 + 𝒲
plot(𝒲, 1e-5)

Aˣ = [1 0; -1 0 ; 0 1.0; 0 -1.0]
bˣ = [0.5; 0.5; 5; 5]
𝕏 = HPolyhedron(Aˣ, bˣ)

# Aᵘ = [1 0; -1 0]
# bᵘ = [15, 15]
# 𝕌 = HPolyhedron(Aᵘ, bᵘ)
𝕌 = Interval(-15.0, 15.0)

# Set tightening for Linear System 2D
function compute_tightened_sets(A, B, K, 𝒲, 𝕏, 𝕌::HPolyhedron; err_approx=5e-6, ϵ=1e-5)
    A_stab = A + B*K  #eigvals(A_stab)
    Ω = approx_mRPI(ϵ, A_stab, 𝒲, err_approx=err_approx)
    𝕏_bar = 𝕏 - Ω
    KΩ = overapproximate(K*Ω, err_approx)
    𝕌_bar = 𝕌 - KΩ
    return 𝕏_bar, 𝕌_bar,  Ω
end

# if 𝕌 is only 1D, i.e. a interval!
function compute_tightened_sets(A, B, K, 𝒲, 𝕏, 𝕌; err_approx=5e-6, ϵ=1e-5)
    A_stab = A + B*K  #eigvals(A_stab)
    Ω = approx_mRPI(ϵ, A_stab, 𝒲, err_approx=err_approx)
    𝕏_bar = 𝕏 - Ω
    KΩ = linear_map(K,Ω)
    𝕌_bar = 𝕌 - KΩ
    !(𝕌_bar ⊆ 𝕌) && error("Interval calculation of 𝕌_bar went wrong!")
    return 𝕏_bar, 𝕌_bar,  Ω
end

@time 𝕏_bar, 𝕌_bar, Ω = compute_tightened_sets(A, B, K, 𝒲, 𝕏, 𝕌)

function l(x1,x2)
    l = x1
    if x2 < -2
        l += 10*x2^2 + 40*x2 + 40
    elseif -2 <= x2 <= 2
        l += 0
    else
        l += 10*x2^2 - 40*x2 + 40
    end
end

## Setup Model

N = 20
n, m = size(B)
z_0 = [-0.258,  3.5]
z_0 = [0,  0]

z_sr = [-0.0076, 0.4275]
z_s = [-0.0358, 2.0009]
z_s = [-0.0675689, 3.78769]
##
mpc_run = Model(with_optimizer(Ipopt.Optimizer))
@variable(mpc_run, z[1:n, 0:N])
@variable(mpc_run, u[1:m, 0:N-1])
constrain_variable_by_set!(mpc_run, z, 𝕏_bar)
constrain_variable_by_set!(mpc_run, u, 𝕌_bar)

# Add dynamics constraints
@constraint(mpc_run, [i=0:N-1, j=1:n],
            z[j, i+1] == A[j,1]*z[1, i] + A[j,2]*z[2, i] + B[j]*u[1,i] )

end_con =     @constraint(mpc_run, [i=1:n], z[i,N] == z_s[i])
initial_con = @constraint(mpc_run, [i=1:n], z[i,0] == z_0[i])

JuMP.register(mpc_run, :l, 2, l, autodiff=true)
@NLobjective(mpc_run, Min, sum(l(z[1,i] - z_s[1], z[2,i] - z_s[2]) for i=0:N) )
# @NLobjective(mpc_run, Min, sum(l(z[1,i], z[2,i]) for i=0:N) )
optimize!(mpc_run)
objective_value(mpc_run)
##

z_real = value.(z).data
u_real = value.(u).data

## Visualize
using Plots
norm_1 = maximum(abs.(z_real[1,:]))
norm_2 = maximum(abs.(z_real[2,:]))
plot(z_real[1,:] , z_real[2,:])
plot!(u_real[1,:])

anim = @gif for i = 1:N
    plot(𝕏_bar, 1e-3, alpha=0.1)
    plot!(z_real[1, :], z_real[2, :], xlim=(-0.5, 0.5), ylim=(-5, 5))
    plot!([z_real[1, i]], [z_real[2, i]], marker=(:hex, 6))
end

list_of_constraint_types(mpc_run)
all_constraints(opt, GenericAffExpr{Float64, VariableRef}, MOI.EqualTo{Float64})
a = all_constraints(mpc_run, GenericAffExpr{Float64, VariableRef}, MOI.LessThan{Float64})

## Animate MPC
z_s = [-0.0675689, 3.78769]
function run_mpc(z_0)
    mpc_run = Model(with_optimizer(Ipopt.Optimizer, print_level=0))
    @variable(mpc_run, z[1:n, 0:N])
    @variable(mpc_run, v[1:m, 0:N-1])
    constrain_variable_by_set!(mpc_run, z, 𝕏_bar)
    constrain_variable_by_set!(mpc_run, v, 𝕌_bar)

    # Add dynamics constraints
    @constraint(mpc_run, [i=0:N-1, j=1:n],
                z[j, i+1] == A[j,1]*z[1, i] + A[j,2]*z[2, i] + B[j]*v[1,i] )

    end_con =     @constraint(mpc_run, [i=1:n], z[i,N] == z_s[i])
    initial_con = @constraint(mpc_run, [i=1:n], z[i,0] == z_0[i])

    JuMP.register(mpc_run, :l, 2, l, autodiff=true)
    @NLobjective(mpc_run, Min, sum(l(z[1,i] - z_s[1], z[2,i] - z_s[2]) for i=0:N) )
    # @NLobjective(mpc_run, Min, sum(l(z[1,i], z[2,i]) for i=0:N) )
    optimize!(mpc_run)
    # objective_value(mpc_run)
    return value.(z).data, value.(v).data
end

# The robot's starting position and velocity
function does()
z= [0,  0]
x = [0,  0]
@gif for i in 1:20
    @show z
    # Plot the current position
    plot(z[1, :], z[2, :], marker=(:hex, 4), xlim=(-0.5, 0.5), ylim=(-5, 5))
    # Run the MPC control optimization
    z_plan, v = run_mpc(z)

    # Draw the planned future states from the MPC optimization
    plot!(z_plan[1, :], z_plan[2, :], linewidth=2)

    # Apply the planned acceleration and simulate one step in time
    v = v[:, 1]
    z⁺ = A*z + B*v
    u = K*(x-z) + v
    x⁺ = A*x + B*u
    @show z⁺, x⁺
    z = z⁺
    x = x⁺
end
end

function simulate_system(A, B, K, 𝒲)

end

does()

plot(𝕏 - Ω + Ω, 1e-3)
plot!(𝕏 - Ω, 1e-3)


text:
lines = readlines("inputs/04.txt") # notes: added a newline to the input file
nums = [parse(Int64, x ) for x in split(lines[1], ',')]
lines = lines[2:end]

function parse_grids(lines)
  grids = []
  for i in 2:6:length(lines)-4
    grid = zeros(Int64, 5,5)
    gridstrs = lines[i:i+4]
    for (r_i, r) in enumerate(gridstrs)
      grid[r_i,:] = [parse(Int64, x) for x in filter(q -> length(q) > 0 ,split(strip(r), ' '))]
    end
    push!(grids, grid)
  end
  return grids
end

is_bingo(grid) = 500 in sum(grid, dims = 1) || 500 in sum(grid, dims = 2)

score(n, grid) = n * (sum(grid) - (100 * count(i -> i == 100, grid)))

function day3(nums, grids_)
  grids = copy(grids_)
  for num in nums
    grids = map(g -> replace(g, num => 100), grids)
    to_del = []
    for (g_i, g) in enumerate(grids)
      if is_bingo(g)
        if length(grids) == 1
          return score(num, grids[1])
        else
          push!(to_del, g_i)
        end
      end
    end
    deleteat!(grids, to_del)
  end
end

grids = parse_grids(lines)
println(day3(nums, grids))


text:
# functions for the parallel particel filter including assiciated help functions

doc"""
    run_pf_paralell(Z::Array{Float64}, theta::Array{Float64}, theta_known::Array{Float64}, N::Int64, N_calc::Float64, dt::Float64, nbr_x0::Int64, nbr_x0_calc::Float64, nbr_x::Int64, nbr_x_calc::Float64, subsample_interval::Int64, subsample_interval_calc::Float64, print_on::Bool, store_weigths::Bool, Κ::Float64, Γ::Float64, A::Float64, B::Float64, c::Float64,d::Float64, f::Float64, g::Float64, power1::Float64, power2::Float64, b_const::Float64)

Runs each of the nbr_pf_proc estiamtions of log-likelihood for the parallel particel filter. Help funciton to pf_paralell.
"""
@fastmath function run_pf_paralell(Z::Array{Float64}, theta::Array{Float64},
  theta_known::Array{Float64}, N::Int64, N_calc::Float64, dt::Float64,dt_U::Float64, nbr_x0::Int64,
  nbr_x0_calc::Float64, nbr_x::Int64, nbr_x_calc::Float64, subsample_interval::Int64,
  subsample_interval_calc::Float64, print_on::Bool, return_weigths_and_particels::Bool,
  Κ::Float64, Γ::Float64, A::Float64, B::Float64, c::Float64,d::Float64, f::Float64,
  g::Float64, power1::Float64, power2::Float64, b_const::Float64)

  # set T
  T = length(Z)

  # Initilized wegits
  w = zeros(N,T)

  # Pre-allocate vector for storing particels
  x = zeros(N,T)

  # set delta for U process
  delta_i = dt_U

  # set start value for nbr_resample
  nbr_resample = 0

  # set start value for loglik
  loglik = 0.

  # pre-allocate vector for logweigths
  logweigths = zeros(N)

  # draw start values for particles
  xinit = 20 + (28-20)*rand(N)

  # pre-allocate vectors
  weigths = zeros(N)
  logweigths = zeros(N)
  dB_x0 = zeros(N,nbr_x0)
  dB_x = zeros(N, nbr_x)
  noise_x = zeros(N,nbr_x)

  # draw start noise
  noise_x0 = randn(N,nbr_x0)

  # set stat value
  xpred = xinit

  # propagate to the first value
  for i = 1:nbr_x0
    @simd for j = 1:N
      #dB_x0[j,i] = sqrt(dt*subsample_interval_calc/nbr_x_calc)*noise_x0[j,i]
      dB_x0[j,i] = sqrt(dt/nbr_x0_calc)*noise_x0[j,i]
    end
  end

  # propagate particels
  propagate_x_v2!(xpred, nbr_x0, nbr_x0_calc, dt, 1, 1., N, N_calc, dB_x0, A, B, c, d, f, g, power1, power2, b_const)

  # store particels
  x[:,1] = xpred

  # calc logweigths
  for j = 1:N
    logweigths[j] =  -log(Γ) - ((Z[1] - x[j,1])/Γ)^2/2 # no unnecessary normalizing constants
  end

  # find max
  constant_w = maximum(logweigths)

  # calc weigths
  for j = 1:N
    weigths[j] = exp(logweigths[j] - constant_w)
  end

  # update loglik
  loglik = loglik + constant_w + log(sum(weigths)) - log(N_calc) # p(Z_0 \mid theta)

  # calc normalized weigths
  w[:,1] = weigths/sum(weigths)

  # update loglik
  #loglik = loglik + log(sum(weigths.*w[:,1])) + sum(constant_w*w[:,1]);

  for t = 2:T

    # fill noise_x with new random numbers
    randn!(noise_x)

    # calc dB_x
    for i = 1:nbr_x
      for j = 1:N
        #dB_x[j,i] = sqrt(dt*subsample_interval_calc/nbr_x_calc)*noise_x[j,i]
        dB_x[j,i] = sqrt(dt/nbr_x_calc)*noise_x[j,i]
      end
    end

    #resample_cond = 1/sum(w[:,t-1].^2) < N/2
    # resample_cond = true # use ordinary boostrap filter
    # resample_cond = 1/sum(w[:,t-1].^2) < N/2  # use a condition for when to resample
    resample_cond = true
    # Only resample if ESS < N/2
    if resample_cond

      # calc resampled indecies
      ind = stratresample(w[:,t-1] , N)

      # resample particels
      x[:,t-1] = x[ind,t-1]

      # update nbr_resample
      nbr_resample = nbr_resample+1

      # reset weigths
      w[:,t-1] = 1/N

    end

    # set xpred
    xpred = x[:,t-1]

    # propagate particles
    propagate_x_v2!(xpred, nbr_x, nbr_x_calc, dt, subsample_interval, subsample_interval_calc, N,N_calc, dB_x,A,B,c,d,f,g,power1,power2,b_const)
    #ccall( (:propagate_x , "./libprop.so"), Void, (Csize_t,Csize_t,Prt{Float64}, Int64,Float64,Float64,Int64,Float64, Int64,Float64, Prt{Float64},Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64),   (N,nbr_x, xpred,nbr_x,nbr_x_calc,dt,subsample_interval,subsample_interval_calc,N,N_calc,dB_x,A,B,c,d,f,g,power1,power2,b_const) )

    # store particels
    x[:,t] = xpred

    # calc logweigths
    for j = 1:N
      logweigths[j] = -log(Γ) - 0.5*log(1-exp(-2*Κ*delta_i))  - (  (Z[t] - xpred[j] - exp(-Κ*delta_i)*(Z[t-1] - x[j,t-1]))/(Γ*sqrt(1-exp(-2*Κ*delta_i)) ) )^2/2.
    end

    # find max
    constant_w = maximum(logweigths)

    # calc weigths
    for j = 1:N
      weigths[j] = w[j,t-1]*exp(logweigths[j] - constant_w) # this is needed since we do not resample every time
    end

    #=
    if ! resample_cond  # use a condition for when to resample

      # if we do not resample we have to "correct" the weigths with w[:,t]
      for j = 1:N
        weigths[j] = w[j,t-1]*weigths[j]
      end

    end
    =#

    # calc normalized weigths
    w[:,t] = weigths/sum(weigths)

    # update loglik
    loglik = loglik + constant_w + log(sum(weigths)) - log(N_calc) # p(Z_t \mid theta)

    #loglik = loglik + log(sum(weigths.*w[:,t])) + sum(constant_w*w[:,t]);

  end

  # calc ESS for the particels at the last time step
  w_end_vec = 1/sum(w[:,end].^2)

  if print_on
    # print ESS and number of resamples
    @printf "ESS: %.4f\n" mean(w_end_vec)
    @printf "Nbr resample: %d \n" round(nbr_resample)
  end

  if return_weigths_and_particels
    return loglik, w, x
  else
    return loglik
  end

end

@inbounds begin

doc"""
    propagate_x_v2!(xpred::Array{Float64}, nbr::Int64,nbr_calc::Float64,dt::Float64,subsample_interval::Int64,subsample_interval_calc::Float64, N::Int64,N_calc::Float64, dB::Array{Float64},A::Float64,B::Float64,c::Float64,d::Float64,f::Float64,g::Float64,power1::Float64,power2::Float64,b_const::Float64)

Propagates the particels one time step.

"""
# @fastmath @inbounds
@fastmath function propagate_x_v2!(xpred::Array{Float64}, nbr::Int64,nbr_calc::Float64,dt::Float64,subsample_interval::Int64,subsample_interval_calc::Float64, N::Int64,N_calc::Float64, dB::Array{Float64},A::Float64,B::Float64,c::Float64,d::Float64,f::Float64,g::Float64,power1::Float64,power2::Float64,b_const::Float64)
  for k = 1:nbr
     @simd for j = 1:N #@simd
      #xpred[j] = xpred[j] - (A*xpred[j] - f + (power2*abs(abs(c - B*xpred[j])^power1/2 - d + g*xpred[j])^(power2 - 1)*sign(abs(c - B*xpred[j])^power1/2 - d + g*xpred[j])*(g - (B*power1*abs(c - B*xpred[j])^(power1 - 1)*sign(c - B*xpred[j]))/2))/2)*dt*subsample_interval_calc/nbr_calc + b_const*dB[j,k]
      xpred[j] = xpred[j] - (A*xpred[j] - f + (power2*abs(abs(c - B*xpred[j])^power1/2 - d + g*xpred[j])^(power2 - 1)*sign(abs(c - B*xpred[j])^power1/2 - d + g*xpred[j])*(g - (B*power1*abs(c - B*xpred[j])^(power1 - 1)*sign(c - B*xpred[j]))/2))/2)*dt/nbr_calc + b_const*dB[j,k]
    end
  end
end

end

doc"""
    stratresample(p , N)

Stratified resampling for sequential Monte Carlo.

Sample N times with repetitions from a distribution on [1:length(p)] with probabilities p. See [link](http://www.cornebise.com/julien/publis/isba2012-slides.pdf).
"""
function stratresample(p , N)

  p = p/sum(p)  # normalize, just in case...

  cdf_approx = cumsum(p)
  cdf_approx[end] = 1
  #I = zeros(N,1)
  indx = zeros(Int64, N)
  U = rand(N,1)
  U = U/N + (0:(N - 1))/N
  index = 1
  for k = 1:N
    while (U[k] > cdf_approx[index])
      index = index + 1
    end
    indx[k] = index
  end

  return indx

end

max_stars_repo_path: DWPSDE model/run_pf_paralell.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: DWPSDE model/run_pf_paralell.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: DWPSDE model/run_pf_paralell.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
include("cpu_stencils.jl")

## Stencil Definition
radius = 4
l = 2*radius + 1
coefs = zeros(Float64, (l,l,l))

for i = 1:radius
    global radius
    posi = radius + i + 1
    negi = radius -i + 1
    mid = radius + 1
    coefs[posi,mid,mid] = round(Float64(1/i),digits=2)
    coefs[negi,mid,mid] = round(Float64(1/i),digits=2)
    coefs[mid, posi, mid] = round(Float64(1/i),digits=2)
    coefs[mid, negi, mid] = round(Float64(1/i),digits=2)
    coefs[mid, mid, posi] = round(Float64(1/i),digits=2)
    coefs[mid, mid, negi] = round(Float64(1/i),digits=2)

end
mid = radius + 1
coefs[mid,mid,mid] = 1
st_def = CreateStencilDefinition(coefs, uses_vsq=true)
st_inst = NewStencilInstance(st_def)
## Input Data size Definition
radius = 4
nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 4
## GPU Implementation
original_data = CreateData(dx,dy,dz)
data = PadData(radius, original_data)
vsq = ones(Float32, dx, dy, dz)
vsq[100:200,:,:] .= 0.1
dout = ApplyStencil(st_inst, data, t_steps,vsq=vsq)
gpu_out = Array(dout)


text:
for f in ("articles", "indefinite_articles", "definite_articles", "prepositions", "pronouns", "stopwords")
  @eval begin
    function $(Symbol(f)){T <: Language}(l::Type{T})
      filename = Pkg.dir("Languages", "data", $f, string(string(l.name.name), ".txt"))
      return fetch_word_list(filename)
    end
  end
end


text:
using Test, RestAPI

const DBFILE = joinpath(dirname(pathof(RestAPI)), "../assets/p2pai.sqlite")
const AUTHFILE = "file://" * joinpath(dirname(pathof(RestAPI)), "../assets/env.dev.json")

# ------------------------------------------------------------ #
# this is demo app, call simple router instead of AUTH_ROUTER
# authentication is skipped
# ------------------------------------------------------------ #
# server = @async RestAPI.run(DBFILE, AUTHFILE)
server = @async RestAPI.init("./assets/p2p.sqlite")

Client.createSysAlert(
    1,
    "SCAL",
    "San Francisco",
    "yellow",
    DateTime("2022-02-23", "yyyy-mm-dd"),
    "ALERT: DocCART Surgical items at par location Oakland are low
in inventory, Auto Replenishment POs must dispatch to avoid deficit."
)

# Client.createUser("amit", "password")
# user = Client.loginUser("amit", "password")

# using HTTP;
# HTTP.CookieRequest.default_cookiejar[1];

# alb1 = Client.createAlbum("Free Yourself Up", "Lake Street Dive", 2018, ["Baby Don't Leave Me Alone With My Thoughts", "Good Kisser"])
# @test Client.pickAlbumToListen() == alb1
# @test Client.pickAlbumToListen() == alb1

# @test Client.getAlbum(alb1.id) == alb1

# push!(alb1.songs, "Shame, Shame, Shame")
# alb2 = Client.updateAlbum(alb1)
# @test length(alb2.songs) == 3
# @test length(Client.getAlbum(alb1.id).songs) == 3

# Client.deleteAlbum(alb1.id)
# Client.pickAlbumToListen()

# alb2 = Client.createAlbum("Haunted Heart", "Charlie Haden Quartet West", 1991, ["Introduction", "Hello My Lovely"])
# @test Client.pickAlbumToListen() == alb2

# Client.createAlbum("Hazards Of Love", "The Decemberists", 2009, ["Prelude", "The Hazards Of Love 1"])
# Client.createAlbum("Tapestry", "Carole King", 1971, ["I Feel The Earth Move", "So Far Away"])


text:
const globalUnique = Ref{Int64}(0)

function genMod(@nospecialize(func), @nospecialize(tt))
    isa(func, Core.Builtin) && error("function is not a generic function")
    world = typemax(UInt)
    meth = which(func, tt)
    sig = Base.signature_type(func, tt)::Type

    (ti, env) = ccall(:jl_type_intersection_with_env, Any,
                      (Any, Any), sig, meth.sig)::Core.SimpleVector
    meth = Base.func_for_method_checked(meth, ti)

    linfo = ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance},
                  (Any, Any, Any, UInt), meth, ti, env, world)

    dependencies = Vector{LLVM.Module}()
    function hook_module_activation(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        push!(dependencies, LLVM.Module(ref))
    end

    params = Base.CodegenParams(cached            = false,
                                prefer_specsig    = true,
                                static_alloc      = false,
                                module_activation = hook_module_activation,
                                )

    # get the code
    mod = let
        ref = ccall(:jl_get_llvmf_defn, LLVM.API.LLVMValueRef,
                    (Any, UInt, Bool, Bool, Base.CodegenParams),
                    linfo, world, #=wrapper=#false, #=optimize=#false, params)
        if ref == C_NULL
            throw(CompilerError(ctx, "the Julia compiler could not generate LLVM IR"))
        end

        llvmf = LLVM.Function(ref)
        LLVM.parent(llvmf)
    end

    return mod, dependencies
end

function irgen(@nospecialize(func), @nospecialize(tt))
    mod, dependencies = genMod(func, tt)

    # the main module should contain a single jfptr_ function definition,
    # e.g. jfptr_kernel_vadd_62977
    definitions = LLVM.Function[]

    for llvmf in functions(mod)
        if !isdeclaration(llvmf)
            push!(definitions, llvmf)
        end
    end

    wrapper = nothing
    for llvmf in definitions
        if startswith(LLVM.name(llvmf), "jfptr_")
            @assert wrapper == nothing
            wrapper = llvmf
        end
    end
    @assert wrapper != nothing

    # the jfptr wrapper function should point us to the actual entry-point,
    # e.g. julia_kernel_vadd_62984
    # FIXME: Julia's globalUnique starting with `-` is probably a bug.
    entry_tag = let
        m = match(r"^jfptr_(.+)_[-\d]+$", LLVM.name(wrapper))
        if m == nothing
            error(LLVM.name(wrapper))
        end
        m.captures[1]
    end
    unsafe_delete!(mod, wrapper)
    entry = let
        re = Regex("^julia_$(entry_tag)_[-\\d]+\$")
        entrypoints = LLVM.Function[]
        for llvmf in definitions
            if llvmf != wrapper
                llvmfn = LLVM.name(llvmf)
                if occursin(re, llvmfn)
                    push!(entrypoints, llvmf)
                end
            end
        end
        if length(entrypoints) != 1
            @warn ":cry:" functions=Tuple(LLVM.name.(definitions)) tag=entry_tag entrypoints=Tuple(LLVM.name.(entrypoints))
        end
        entrypoints[1]
    end

    # link in dependent modules
    for dep in dependencies
        link!(mod, dep)
    end

    # remove other fptr wrappers
    for llvmf in functions(mod)
        if !isdeclaration(llvmf) && startswith(LLVM.name(llvmf), "jfptr_")
            unsafe_delete!(mod, llvmf)
        end
    end

    # rename the entry point
    llvmfn = replace(LLVM.name(entry), r"_\d+$"=>"")

    ## append a global unique counter
    globalUnique[] += 1
    llvmfn *= "_$(globalUnique[])"
    LLVM.name!(entry, llvmfn)

    return mod, entry
end


text:
#
# sparse-benchs.jl -
#
# Benchmarking sparse operators.
#
module BenchmarkingLazyAlgebraSparseOperators

using LinearAlgebra
using SparseArrays
using StructuredArrays
using LazyAlgebra
using LazyAlgebra: identical, Adjoint, Direct
using LazyAlgebra.SparseMethods
using LazyAlgebra.SparseOperators: check_structure, compute_offsets
using BenchmarkTools
using Random

is_csc(::Any) = false
is_csc(::CompressedSparseOperator{:CSC}) = true
is_csc(::Adjoint{<:CompressedSparseOperator{:CSR}}) = true

is_csr(::Any) = false
is_csr(::CompressedSparseOperator{:CSR}) = true
is_csr(::Adjoint{<:CompressedSparseOperator{:CSC}}) = true

is_coo(::Any) = false
is_coo(::CompressedSparseOperator{:COO}) = true
is_coo(::Adjoint{<:CompressedSparseOperator{:COO}}) = true

# Generate a possibly sparse array of random values.  Value are small signed
# integers so that all computations should be exact (except with non-integer
# multipliers).
genarr(T::Type, dims::Integer...; kwds...) = genarr(T, dims; kwds...)
function genarr(T::Type, dims::Tuple{Vararg{Integer}};
                sparsity::Real = 0,
                range::AbstractUnitRange{<:Integer} = -17:17)
    @assert 0 ≤ sparsity ≤ 1
    A = Array{T}(undef, dims)
    for i in eachindex(A)
        if sparsity > 0 && rand() ≤ sparsity
            A[i] = zero(T)
        elseif T <: Complex
            A[i] = T(rand(range), rand(range))
        else
            A[i] = T(rand(range))
        end
    end
    return A
end

# Unpack a sparse operator into a regular array using simplest iterator.  There
# may be duplicates.
function unpack_with_iterator!(dest::Array{T},
                               A::SparseOperator{E},
                               op = (E === Bool ? (|) : (+))) where {T,E}
    C = fill!(reshape(dest, (nrows(A), ncols(A))), zero(T))
    for (Aij, i, j) in A
        C[i,j] = op(C[i,j], Aij)
    end
    return dest
end

nickname(::AbstractMatrix) = "Matrix";
nickname(::AbstractSparseMatrix) = "SparseMatrix";
nickname(::GeneralMatrix) = "GeneralMatrix";
nickname(::CompressedSparseOperator{:COO}) = "SparseOperatorCOO";
nickname(::CompressedSparseOperator{:CSC}) = "SparseOperatorCSC";
nickname(::CompressedSparseOperator{:CSR}) = "SparseOperatorCSR";
nickname(n::Integer, x) = rpad(nickname(x), n);

function bench1(; m::Integer = 2_000, n::Integer = m, T::Type = Float64,
                sparsity::Real = 0.95)
    bench1(T, m, n, sparsity)
end

bench1(::Type{T}, m::Integer, n::Integer, sparsity::Real) where {T} =
    bench1(T, Int(m), Int(n), Float64(sparsity))

function bench1(::Type{T}, m::Int, n::Int, sparsity::Float64) where {T}
    A = genarr(T, (m, n); sparsity=sparsity);
    x = genarr(T, n);
    y = genarr(T, m);
    x1 = copy(x);
    y0 = copy(y);
    x1 = similar(x);
    y1 = similar(y);
    x2 = similar(x);
    y2 = similar(y);
    coo = SparseOperatorCOO(A);
    csc = SparseOperatorCSC(A);
    csr = SparseOperatorCSR(A);
    gen = GeneralMatrix(A);
    S = sparse(A);
    println("Tests are done for T=$T, (m,n)=($m,$n) and sparsity = ",
            round(sparsity*1e2, sigdigits=3), "% of entries.\n")
    nnz(S) == nnz(coo) || println("not same number on non-zeros (COO)");
    nnz(S) == nnz(csc) || println("not same number on non-zeros (CSS)");
    nnz(S) == nnz(csr) || println("not same number on non-zeros (CSR)");
    mul!(y1, A, x);
    mul!(x1, A', y);
    for B in (gen, S, coo, csc, csr)
        mul!(y2, B, x);
        mul!(x2, B', y);
        println("compare A*x  for a ", nickname(20, B), extrema(x1 - x2))
        println("compare A'*x for a ", nickname(20, B), extrema(y1 - y2))
    end
    println()
    for B in (A, gen, S, coo, csc, csr)
        print("benchmarking A*x for a ", nickname(20, B))
        @btime mul!($y1, $B, $x)
    end
    println()
    for B in (A, gen, S, coo, csc, csr)
        print("benchmarking A'*x for a ", nickname(20, B))
        @btime mul!($x1, $(B'), $y)
    end

end

end # module

BenchmarkingLazyAlgebraSparseOperators.bench1()


text:
####################
# PullRequest Type #
####################

mutable struct PullRequest <: GitHubType
    base::Union{Branch, Nothing}
    head::Union{Branch, Nothing}
    number::Union{Int, Nothing}
    id::Union{Int, Nothing}
    comments::Union{Int, Nothing}
    commits::Union{Int, Nothing}
    additions::Union{Int, Nothing}
    deletions::Union{Int, Nothing}
    changed_files::Union{Int, Nothing}
    state::Union{String, Nothing}
    title::Union{String, Nothing}
    body::Union{String, Nothing}
    merge_commit_sha::Union{String, Nothing}
    created_at::Union{Dates.DateTime, Nothing}
    updated_at::Union{Dates.DateTime, Nothing}
    closed_at::Union{Dates.DateTime, Nothing}
    merged_at::Union{Dates.DateTime, Nothing}
    url::Union{HTTP.URI, Nothing}
    html_url::Union{HTTP.URI, Nothing}
    assignee::Union{Owner, Nothing}
    user::Union{Owner, Nothing}
    merged_by::Union{Owner, Nothing}
    milestone::Union{Dict, Nothing}
    _links::Union{Dict, Nothing}
    mergeable::Union{Bool, Nothing}
    merged::Union{Bool, Nothing}
    locked::Union{Bool, Nothing}
end

PullRequest(data::Dict) = json2github(PullRequest, data)
PullRequest(number::Real) = PullRequest(Dict("number" => number))
namefield(pr::PullRequest) = pr.number

mutable struct PullRequestFile <: GitHubType
    raw_url::Union{String, Nothing}
    status::Union{String, Nothing}
    patch::Union{String, Nothing}
    changes::Union{Int, Nothing}
    sha::Union{String, Nothing}
    filename::Union{String, Nothing}
    additions::Union{Int, Nothing}
    deletions::Union{Int, Nothing}
    blob_url::Union{String, Nothing}
    contents_url::Union{String, Nothing}
end

PullRequestFile(data::Dict) = json2github(PullRequestFile, data)
PullRequestFile(fname::String) = PullRequestFile(Dict("filename" => fname))
namefield(prf::PullRequestFile) = prf.filename

###############
# API Methods #
###############

@api_default function pull_requests(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/pulls"; options...)
    return map(PullRequest, results), page_data
end

@api_default function pull_request(api::GitHubAPI, repo, pr; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/pulls/$(name(pr))"; options...)
    return PullRequest(result)
end

@api_default function pull_request_files(api::GitHubAPI, repo, pr; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/pulls/$(name(pr))/files"; options...)
    return [PullRequestFile(f) for f in result]
end

@api_default function update_pull_request(api::GitHubAPI, repo, pr; options...)
    result = gh_patch_json(api, "/repos/$(name(repo))/pulls/$(name(pr))"; options...)
    return PullRequest(result)
end

@api_default function close_pull_request(api::GitHubAPI, repo, pr; options...)
    update_pull_request(api, repo, pr; params = Dict(
        :state => "closed"
    ), options...)
end

@api_default function merge_pull_request(api::GitHubAPI, repo, pr; options...)
    gh_put_json(api, "/repos/$(name(repo))/pulls/$(name(pr))/merge"; options...)
end

@api_default function create_pull_request(api::GitHubAPI, repo; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/pulls"; options...)
    return PullRequest(result)
end

@api_default function create_comment(api::GitHubAPI, repo, pr::PullRequest, body::AbstractString; options...)
    create_comment(api, repo, pr, :pr; params = Dict(
        :body => body
    ), options...)
end


text:
"""
    $(TYPEDEF)

A struct used to reference the `equilibrator_api` object.

$(FIELDS)
"""
struct Equilibrator
    eq::PyObject
    cc::PyObject
end

function Base.show(io::IO, ::MIME"text/plain", equilibrator::Equilibrator)
    v = string(equilibrator.eq.__version__)
    println(crayon"blue bold", "Interface to equilibrator_api v$v", crayon"default !bold")
    println(
        crayon"yellow",
        "Temperature:\t",
        crayon"default",
        get_temperature(equilibrator),
    )
    println(
        crayon"yellow",
        "Ionic strength: ",
        crayon"default",
        get_ionic_strength(equilibrator),
    )
    println(crayon"yellow", "pH:\t\t", crayon"default", get_pH(equilibrator))
    println(crayon"yellow", "pMg:\t\t", crayon"default", get_pMg(equilibrator))
end

"""
    $(SIGNATURES)

A constructor for `Equilibrator`. Keyword arguments `ionic_strength` and `temperature` have units.
Any appropriate unit may be used, the conversion is handled internally. This is typically
slow to load, ~ 30 seconds, since all the thermodynamic data needs to be loaded into memory.

# Example
```
using eQuilibrator
using Unitful

equilibrator = Equilibrator(pH=6.0, pMg = 2.5, ionic_strength=250.0u"mM", temperature=30.0u"°C")
```
"""
function Equilibrator(;
    pH = 7.5,
    pMg = 3.0,
    ionic_strength = 0.25u"M",
    temperature = 298.15u"K",
)
    eq = pyimport("equilibrator_api")
    cc = eq.ComponentContribution()

    cc.p_h = eq.Q_(pH)
    cc.p_mg = eq.Q_(pMg)

    _ionic_strength = uconvert(u"M", float(ionic_strength))
    cc.ionic_strength = eq.Q_(string(_ionic_strength))

    _temperature = uconvert(u"K", float(temperature))
    cc.temperature = eq.Q_(string(_temperature))
    return Equilibrator(eq, cc)
end

"""
    $(SIGNATURES)

Reset `equilibrator` back to the default starting state with `pH = 7.5`,
`pMg = 3.0`, `temperature = 298 K` and `ionic_strength = 0.25 M`.
"""
function reset_equilibrator(equilibrator::Equilibrator)
    set_temperature(equilibrator, 298.15u"K")
    set_pH(equilibrator, 7.5)
    set_pMg(equilibrator, 3.0)
    set_ionic_strength(equilibrator, 0.25u"M")
end

"""
    $(SIGNATURES)

Get the temperature of `equilibrator`.

# Example
```
get_temperature(equilibrator)
```
"""
function get_temperature(equilibrator::Equilibrator)
    units = uparse(_parse_units(py"str"(equilibrator.cc.temperature.units)))
    (equilibrator.cc.temperature.m)units
end

"""
    $(SIGNATURES)

Set the temperature of `equilibrator` to `x`.

# Example
```
set_temperature(equilibrator, 40u"°C")
```
"""
function set_temperature(equilibrator::Equilibrator, x)
    _temperature = uconvert(u"K", float(x))
    equilibrator.cc.temperature = equilibrator.eq.Q_(string(_temperature))
    return nothing
end

"""
    $(SIGNATURES)

Get the ionic strength of `equilibrator`.

# Example
```
get_ionic_strength(equilibrator)
```
"""
function get_ionic_strength(equilibrator::Equilibrator)
    units = uparse(_parse_units(py"str"(equilibrator.cc.ionic_strength.units)))
    (equilibrator.cc.ionic_strength.m)units
end

"""
    $(SIGNATURES)

Set the ionic strength of `equilibrator` to `x`.

# Example
```
set_ionic_strength(equilibrator, 100u"mM")
```
"""
function set_ionic_strength(equilibrator::Equilibrator, x)
    _ionic_strength = uconvert(u"M", float(x))
    equilibrator.cc.ionic_strength = equilibrator.eq.Q_(string(_ionic_strength))
    return nothing
end

"""
    $(SIGNATURES)

Get the pH from `equilibrator`.

# Example
```
get_pH(equilibrator)
```
"""
function get_pH(equilibrator::Equilibrator)
    equilibrator.cc.p_h.m
end

"""
    $(SIGNATURES)

Set pH of `equilibrator` to `x`.

# Example
```
set_pH(equilibrator, 3.5)
```
"""
function set_pH(equilibrator::Equilibrator, x)
    equilibrator.cc.p_h = equilibrator.eq.Q_(float(x))
    return nothing
end

"""
    $(SIGNATURES)

Get pMg from `equilibrator`.

# Example
```
get_pMg(equilibrator)
```
"""
function get_pMg(equilibrator::Equilibrator)
    equilibrator.cc.p_mg.m
end

"""
    $(SIGNATURES)

Set pMg of `equilibrator` to `x`.

# Example
```
set_pMg(equilibrator, 1.2)
```
"""
function set_pMg(equilibrator::Equilibrator, x)
    equilibrator.cc.p_mg = equilibrator.eq.Q_(float(x))
    return nothing
end


text:
function resolve(path)
    if islink(path)
        dir = dirname(path)
        resolve(joinpath(dir, readlink(path)))
    else
        path
    end
end

# return a list of valid directories, resolving symlinks and pruning duplicates
function valid_dirs(dirs)
    map!(resolve, dirs, dirs)
    filter(isdir, unique(dirs))
end

function join_versions(versions)
    isempty(versions) && return "no specific version"
    "version " * join(versions, " or ")
end

function join_locations(locations)
    isempty(locations) && return "in no specific location"
    "in " * join(locations, " or ")
end

## generic discovery routines

function library_names(name::String, versions::Vector=[])
    names = String[]

    # always look for an unversioned library first
    if Sys.iswindows()
        push!(names, "$(name)$(Sys.WORD_SIZE).$(Libdl.dlext)")

        # some libraries (e.g. CUTENSOR) are shipped without the word size-prefix
        push!(names, "$(name).$(Libdl.dlext)")
    elseif Sys.isapple()
        # macOS puts the version number before the dylib extension
        push!(names, "lib$(name).$(Libdl.dlext)")
    elseif Sys.isunix()
        # most UNIX distributions ship versioned libraries (also see JuliaLang/julia#22828)
        push!(names, "lib$(name).$(Libdl.dlext)")
    else
        push!(names, "lib$name.$(Libdl.dlext)")
    end

    # then consider versioned libraries
    for version in versions
        if Sys.iswindows()
            # Windows encodes the version in the filename
            if version isa VersionNumber
                append!(names, ["$(name)$(Sys.WORD_SIZE)_$(version.major)$(version.minor).$(Libdl.dlext)",
                                "$(name)$(Sys.WORD_SIZE)_$(version.major).$(Libdl.dlext)"])
            elseif version isa String
                push!(names, "$(name)$(Sys.WORD_SIZE)_$(version).$(Libdl.dlext)")
            end

            # some libraries (e.g. CUTENSOR) are shipped without the word size-prefix
            if version isa VersionNumber
                append!(names, ["$(name)_$(version.major)$(version.minor).$(Libdl.dlext)",
                                "$(name)_$(version.major).$(Libdl.dlext)"])
            elseif version isa String
                push!(names, "$(name)_$(version).$(Libdl.dlext)")
            end
        elseif Sys.isapple()
            # macOS puts the version number before the dylib extension
            if version isa VersionNumber
                append!(names, ["lib$(name).$(version.major).$(version.minor).$(Libdl.dlext)",
                                "lib$(name).$(version.major).$(Libdl.dlext)"])
            elseif version isa String
                push!(names, "lib$(name).$(version).$(Libdl.dlext)")
            end
        elseif Sys.isunix()
            # most UNIX distributions ship versioned libraries (also see JuliaLang/julia#22828)
            if version isa VersionNumber
                append!(names, ["lib$(name).$(Libdl.dlext).$(version.major).$(version.minor).$(version.patch)",
                                "lib$(name).$(Libdl.dlext).$(version.major).$(version.minor)",
                                "lib$(name).$(Libdl.dlext).$(version.major)"])
            elseif version isa String
                push!(names, "lib$(name).$(Libdl.dlext).$(version)")
            end
        end
    end

    return names
end

"""
    find_library(name, versions; locations=String[])

Wrapper for Libdl.find_library, performing a more exhaustive search:

- variants of the library name (including version numbers, platform-specific tags, etc);
- various subdirectories of the `locations` list, and finally system library directories.

Returns the full path to the library.
"""
function find_library(name::String, versions::Vector=[];
                      locations::Vector{String}=String[])
    # figure out names
    all_names = library_names(name, versions)

    # figure out locations
    all_locations = String[]
    for location in locations
        push!(all_locations, location)
        push!(all_locations, joinpath(location, "lib"))
        if Sys.WORD_SIZE == 64
            push!(all_locations, joinpath(location, "lib64"))
            push!(all_locations, joinpath(location, "libx64"))
        end
        if Sys.iswindows()
            push!(all_locations, joinpath(location, "bin"))
            push!(all_locations, joinpath(location, "bin", Sys.WORD_SIZE==64 ? "x64" : "Win32"))
        end
    end

    @debug "Looking for library $name, $(join_versions(versions)), $(join_locations(locations))" all_names all_locations
    name_found = Libdl.find_library(all_names, all_locations)
    if isempty(name_found)
        @debug "Did not find $name"
        return nothing
    end

    # find the full path of the library (which Libdl.find_library doesn't guarantee to return)
    path = Libdl.dlpath(name_found)
    @debug "Found $(basename(path)) at $(dirname(path))"
    return path
end

"""
    find_binary(name; locations=String[])

Similar to `find_library`, performs an exhaustive search for a binary in various
subdirectories of `locations`, and finally PATH by using `Sys.which`.
"""
function find_binary(name::String; locations::Vector{String}=String[])
    # figure out locations
    all_locations = String[]
    for location in locations
        push!(all_locations, location)
        push!(all_locations, joinpath(location, "bin"))
    end
    # we look in PATH too by using `Sys.which` with unadorned names

    @debug "Looking for binary $name $(join_locations(locations))" all_locations
    all_paths = [name; [joinpath(location, name) for location in all_locations]]
    for path in all_paths
        try
            program_path = Sys.which(path)
            if program_path !== nothing
                @debug "Found $path at $program_path"
                return program_path
            end
        catch
            # some system disallow `stat` on certain paths
        end
    end

    @debug "Did not find $path"
    return nothing
end

## CUDA-specific discovery routines

const cuda_releases = [v"9.0", v"9.1", v"9.2",
                       v"10.0", v"10.1", v"10.2",
                       v"11.0", v"11.1", v"11.2", v"11.3", v"11.4"]

# return possible versions of a CUDA library
function cuda_library_versions(name::String)
    if Sys.iswindows()
        # CUDA libraries on Windows are always versioned, however, we don't
        # know which version we're looking for (and we don't first want to
        # figure that out by, say, invoking a versionless binary like ptxas).

        # start out with all known CUDA releases
        versions = Any[cuda_releases...]

        # append some future releases
        for major in last(versions).major:15, minor in 1:10
            version = VersionNumber(major, minor)
            if !in(version, versions)
                push!(versions, version)
            end
        end

        # CUPTI is special, and uses a dot-separated, year-based versioning
        if name == "cupti"
            for year in 2020:2022, major in 1:5, minor in 0:3
                version = "$year.$major.$minor"
                push!(versions, version)
            end
        end

        # NVTX is special, and only uses a single digit
        if name == "nvToolsExt"
            append!(versions, [v"1", v"2"])
        end

        versions
    else
        # only consider unversioned libraries on other platforms.
        []
    end
end

# simplified find_library/find_binary entry-points,
# looking up name aliases and known version numbers
# and passing the (optional) toolkit dirs as locations.
function find_cuda_library(toolkit_dirs::Vector{String}, library::String, versions::Vector)
    # figure out the location
    locations = toolkit_dirs
    ## CUPTI is in the "extras" directory of the toolkit
    if library == "cupti"
        toolkit_extras_dirs = filter(dir->isdir(joinpath(dir, "extras")), toolkit_dirs)
        cupti_dirs = map(dir->joinpath(dir, "extras", "CUPTI"), toolkit_extras_dirs)
        append!(locations, cupti_dirs)
    end
    ## NVTX is located in an entirely different location on Windows
    if library == "nvToolsExt" && Sys.iswindows()
        if haskey(ENV, "NVTOOLSEXT_PATH")
            dir = ENV["NVTOOLSEXT_PATH"]
            @debug "Looking for NVTX library via environment variable" dir
        else
            program_files = ENV[Sys.WORD_SIZE == 64 ? "ProgramFiles" : "ProgramFiles(x86)"]
            dir = joinpath(program_files, "NVIDIA Corporation", "NvToolsExt")
            @debug "Looking for NVTX library in the default directory" dir
        end
        isdir(dir) && push!(locations, dir)
    end

    find_library(library, versions; locations)
end
function find_cuda_binary(toolkit_dirs::Vector{String}, name::String)
    # figure out the location
    locations = toolkit_dirs
    ## compute-sanitizer is in the "extras" directory of the toolkit
    if name == "compute-sanitizer"
        toolkit_extras_dirs = filter(dir->isdir(joinpath(dir, "extras")), toolkit_dirs)
        sanitizer_dirs = map(dir->joinpath(dir, "extras", "compute-sanitizer"), toolkit_extras_dirs)
        append!(locations, sanitizer_dirs)
    end

    find_binary(name; locations)
end

"""
    find_toolkit()::Vector{String}

Look for directories where (parts of) the CUDA toolkit might be installed. This returns a
(possibly empty) list of paths that can be used as an argument to other discovery functions.

The behavior of this function can be overridden by defining the `CUDA_PATH`, `CUDA_HOME` or
`CUDA_ROOT` environment variables, which should point to the root of the CUDA toolkit.
"""
function find_toolkit()
    dirs = String[]

    # look for environment variables to override discovery
    envvars = ["CUDA_PATH", "CUDA_HOME", "CUDA_ROOT"]
    filter!(var -> haskey(ENV, var) && ispath(ENV[var]), envvars)
    if !isempty(envvars)
        paths = unique(map(var->ENV[var], envvars))
        if length(paths) > 1
            @warn "Multiple CUDA environment variables set to different values: $(join(paths, ", "))"
        end

        @debug "Looking for CUDA toolkit via environment variables $(join(envvars, ", "))"
        append!(dirs, paths)
        return dirs
    end

    # look for the compiler binary (in the case PATH points to the installation)
    ptxas_path = find_binary("ptxas")
    if ptxas_path !== nothing
        ptxas_dir = dirname(ptxas_path)
        if occursin(r"^bin(32|64)?$", basename(ptxas_dir))
            ptxas_dir = dirname(ptxas_dir)
        end

        @debug "Looking for CUDA toolkit via ptxas binary" path=ptxas_path dir=ptxas_dir
        push!(dirs, ptxas_dir)
    end

    # look for the runtime library (in the case LD_LIBRARY_PATH points to the installation)
    libcudart_path = find_library("cudart")
    if libcudart_path !== nothing
        libcudart_dir = dirname(libcudart_path)
        if occursin(r"^(lib|bin)(32|64)?$", basename(libcudart_dir))
            libcudart_dir = dirname(libcudart_dir)
        end

        @debug "Looking for CUDA toolkit via CUDA runtime library" path=libcudart_path dir=libcudart_dir
        push!(dirs, libcudart_dir)
    end

    # look in default installation directories
    default_dirs = String[]
    if Sys.iswindows()
        # CUDA versions are installed in separate directories under a single base dir
        program_files = ENV[Sys.WORD_SIZE == 64 ? "ProgramFiles" : "ProgramFiles(x86)"]
        basedir = joinpath(program_files, "NVIDIA GPU Computing Toolkit", "CUDA")
        if isdir(basedir)
            entries = map(dir -> joinpath(basedir, dir), readdir(basedir))
            append!(default_dirs, entries)
        end
    else
        # CUDA versions are installed in unversioned dirs, or suffixed with the version
        basedirs = ["/usr/local/cuda", "/opt/cuda"]
        for ver in cuda_releases, dir in basedirs
            push!(default_dirs, "$dir-$(ver.major).$(ver.minor)")
        end
        append!(default_dirs, basedirs)
        push!(default_dirs, "/usr/lib/nvidia-cuda-toolkit")
        push!(default_dirs, "/usr/share/cuda")
    end
    reverse!(default_dirs) # we want to search starting from the newest CUDA version
    default_dirs = valid_dirs(default_dirs)
    if !isempty(default_dirs)
        @debug "Looking for CUDA toolkit via default installation directories" dirs=default_dirs
        append!(dirs, default_dirs)
    end

    # filter
    dirs = valid_dirs(dirs)
    @debug "Found CUDA toolkit at $(join(dirs, ", "))"
    return dirs
end

"""
    find_libdevice(toolkit_dirs::Vector{String})

Look for the CUDA device library supporting `targets` in any of the CUDA toolkit directories
`toolkit_dirs`. On CUDA >= 9.0, a single library unified library is discovered and returned
as a string. On older toolkits, individual libraries for each of the targets are returned as
a vector of strings.
"""
function find_libdevice(toolkit_dirs)
    @debug "Request to look for libdevice" locations=toolkit_dirs

    # figure out locations
    dirs = String[]
    for toolkit_dir in toolkit_dirs
        push!(dirs, toolkit_dir)
        push!(dirs, joinpath(toolkit_dir, "libdevice"))
        push!(dirs, joinpath(toolkit_dir, "nvvm", "libdevice"))
    end

    # filter
    dirs = valid_dirs(dirs)
    @debug "Look for libdevice" locations=dirs

    for dir in dirs
        path = joinpath(dir, "libdevice.10.bc")
        if isfile(path)
            @debug "Found unified device library at $path"
            return path
        end
    end

    return nothing
end

"""
    find_libcudadevrt(toolkit_dirs::Vector{String})

Look for the CUDA device runtime library in any of the CUDA toolkit directories
`toolkit_dirs`.
"""
function find_libcudadevrt(toolkit_dirs)
    locations = toolkit_dirs
    @debug "Request to look for libcudadevrt " locations

    name = nothing
    if Sys.isunix()
        name = "libcudadevrt.a"
    elseif Sys.iswindows()
        name = "cudadevrt.lib"
    else
        error("No support for discovering the CUDA device runtime library on your platform, please file an issue.")
    end

    # figure out locations
    all_locations = String[]
    for location in locations
        push!(all_locations, location)
        if Sys.iswindows()
            if Sys.WORD_SIZE == 64
                push!(all_locations, joinpath(location, "lib", "x64"))
            elseif Sys.WORD_SIZE == 32
                push!(all_locations, joinpath(location, "lib", "Win32"))
            end
        else
            push!(all_locations, joinpath(location, "lib"))
            if Sys.WORD_SIZE == 64
                push!(all_locations, joinpath(location, "lib64"))
            end
        end
    end

    @debug "Looking for CUDA device runtime library $name" locations=all_locations
    paths = filter(isfile, map(location->joinpath(location, name), all_locations))

    if isempty(paths)
        return nothing
    else
        path = first(paths)
        @debug "Found CUDA device runtime library $(basename(path)) at $(dirname(path))"
        return path
    end
end