text:
module plate_w_hole_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshExportModule
using LinearAlgebra: norm
using SparseArrays: cholesky
using Statistics: mean

function plate_w_hole_H20_stress()
    E = 2.4*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    Re = 0.3*phun("M"); # outer radius
    Ri= 0.1*phun("M"); # hole radius
    H = 0.1*phun("M") # thickness of the plate
    nRadial, nCircumferential, nThickness = 6, 8, 1;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    convergencestudy = FDataDict[]
    for ref in 0:1:2
        println("ref = $(ref)")
        # Thickness = H
        Thickness = H/2^ref
        tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

        fens,fes = H20block(1.0, pi/2, Thickness, 2^ref*nRadial, 2^ref*nCircumferential, 2^ref*nThickness)

        bdryfes = meshboundary(fes);
        icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);

        for i=1:count(fens)
            t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
            fens.xyz[i,:] = [(t*Re+(1-t)*Ri)*cos(a), (t*Re+(1-t)*Ri)*sin(a), z];
        end

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

        l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
        setebc!(u,l1,true, 2, 0.0)
        l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
        setebc!(u,l1,true, 1, 0.0)
        # Plane-stress constraint: assume the plane z=0 is the plane of symmetry of the plate
        l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
        setebc!(u,l1,true, 3, 0.0)
        # If this was enabled, the plane-strain  constraint would be enforced.
        # l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
        # setebc!(u,l1,true, 3, 0.0)

        applyebc!(u)
        numberdofs!(u)
        el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), GaussRule(2, 3)))
        function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
            local r = sqrt(XYZ[1]^2 + XYZ[2]^2)
            nx = XYZ[1]/r; ny = XYZ[2]/r
            forceout[1] = sigmarr(XYZ) * nx - sigmart(XYZ) * ny
            forceout[2] = sigmarr(XYZ) * ny + sigmart(XYZ) * nx
            forceout[3] = 0.0
            return forceout
        end
        fi = ForceIntensity(FFlt, 3, pfun);
        F2 = distribloads(el1femm, geom, u, fi, 2);

        MR = DeforModelRed3D

        material = MatDeforElastIso(MR, E, nu)

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

        # The geometry field now needs to be associated with the FEMM
        femm = associategeometry!(femm, geom)

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

        stressfields = elemfieldfromintegpoints(femm, geom, u, :Cauchy, collect(1:6))

        push!(convergencestudy, FDataDict(
        "elementsize"=> 1.0 / 2^ref,
        "fens"=>fens,
        "fes"=>fes,
        "geom"=>geom,
        "u"=>u,
        "femm"=>femm,
        "stressfields"=>[stressfields],
        "tolerance"=>tolerance)
        )
    end # for ref in

    # File = "mplate_w_hole_H20m_stress"
    # open(File * ".jls", "w") do file
    #     serialize(file, convergencestudy)
    # end

end # plate_w_hole_H20_stress

function plate_w_hole_MSH8_convergence()
    E = 2.4*phun("MEGA*PA");# 210e3 MPa
    nu = 0.49995;
    Re = 0.3*phun("M"); # outer radius
    Ri= 0.1*phun("M"); # hole radius
    H = 0.1*phun("M") # thickness of the plate
    nRadial, nCircumferential=3, 5;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    sigxderrs = Dict{Symbol, FFltVec}()
    sigyderrs = Dict{Symbol, FFltVec}()
    numelements = []
    numnodes = []
    for extrap in [:extraptrend :extrapmean]
        sigxderrs[extrap] = FFltVec[]
        sigyderrs[extrap] = FFltVec[]
        numelements = []
        numnodes = []
        for ref in 0:1:2
            # Thickness = H
            Thickness = H/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = H8block(1.0, pi/2, Thickness, 2^ref*nRadial, 2^ref*nCircumferential, 1)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);

            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
                fens.xyz[i,:] = [(t*Re+(1-t)*Ri)*cos(a), (t*Re+(1-t)*Ri)*sin(a), z];
            end

            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, fens.xyz,
            #     FinEtools.MeshExportModule.H8)
            # @async run(`"paraview.exe" $File`)

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

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            # Plane-stress constraint: assume the plane z=0 is the plane of symmetry of the plate
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)
            # If this was enabled, the plane-strain  constraint would be enforced.
            # l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
            # setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)
            el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), GaussRule(2, 2)))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                local r = sqrt(XYZ[1]^2 + XYZ[2]^2)
                nx = XYZ[1]/r; ny = XYZ[2]/r
                # local sx, sy, txy
                # sx, sy, txy = sigmaxx(XYZ), sigmayy(XYZ), sigmaxy(XYZ)
                # sn = sx * nx^2 + sy * ny^2 + 2 * nx * ny * txy
                # tn = -(sx - sy) * nx * ny + (nx^2 - ny^2) * txy
                # forceout[1] = sn * nx - tn * ny
                # forceout[2] = sn * ny + tn * nx
                # forceout[3] = 0.0
                forceout[1] = sigmarr(XYZ) * nx - sigmart(XYZ) * ny
                forceout[2] = sigmarr(XYZ) * ny + sigmart(XYZ) * nx
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

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

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

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

            nlA = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            nlB = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, Thickness], inflate=tolerance);
            # thecorneru = zeros(FFlt,length(nlA),3)
            # gathervalues_asmat!(u, thecorneru, nl);
            # thecorneru = mean(thecorneru, 1)[1]/phun("mm")
            # println("displacement = $(thecorneru) vs -0.10215 [MM]")

            println("Extrapolation: $( extrap )--------------- ")
            sigx = fieldfromintegpoints(femm, geom, u, :Cauchy, 1;
            nodevalmethod = :averaging, reportat = extrap)
            sigy = fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
            nodevalmethod = :averaging, reportat = extrap)
            sigyA = mean(sigy.values[nlA,1], dims = 1)[1]
            sigyAtrue = sigmatt([Ri, 0.0, 0.0])
            println("sig_y@A =$(sigyA/phun("MPa")) vs $(sigyAtrue/phun("MPa")) [MPa]")
            sigxB = mean(sigx.values[nlB,1], dims = 1)[1]
            sigxBtrue = sigmatt([0.0, Ri, 0.0])
            println("sig_x@B =$(sigxB/phun("MPa")) vs $(sigxBtrue/phun("MPa")) [MPa]")
            push!(numnodes, count(fens))
            push!(numelements, count(fes))
            push!(sigxderrs[extrap], abs(sigxB/sigxBtrue - 1.0))
            push!(sigyderrs[extrap], abs(sigyA/sigyAtrue - 1.0))
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            # FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
            # scalars=[("sigmax", sigx.values/phun("MEGA*PA"))])
            # @async run(`"paraview.exe" $File`)
        end
    end

    File = "plate_w_hole_MSH8_convergence.CSV"
    savecsv(File, numelements=vec(numelements), numnodes=vec(numnodes), sigxderrtrend=vec(sigxderrs[:extraptrend]), sigxderrdefault=vec(sigxderrs[:extrapmean]), sigyderrtrend=vec(sigyderrs[:extraptrend]), sigyderrdefault=vec(sigyderrs[:extrapmean]))

end # plate_w_hole_MSH8_convergence

function plate_w_hole_MSH8_PE_convergence()
    E = 2.4*phun("MEGA*PA");# 210e3 MPa
    nu = 0.49995;
    Re = 0.3*phun("M"); # outer radius
    Ri= 0.1*phun("M"); # hole radius
    H = 0.1*phun("M") # thickness of the plate
    nRadial, nCircumferential=3, 5;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    sigxderrs = Dict{Symbol, FFltVec}()
    sigyderrs = Dict{Symbol, FFltVec}()
    numelements = []
    numnodes = []
    for extrapolation in [:extraptrend :extrapmean]
        sigxderrs[extrapolation] = FFltVec[]
        sigyderrs[extrapolation] = FFltVec[]
        numelements = []
        numnodes = []
        for ref in 0:1:2
            # Thickness = H
            Thickness = H/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = H8block(1.0, pi/2, Thickness, 2^ref*nRadial, 2^ref*nCircumferential, 1)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);

            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
                fens.xyz[i,:] = [(t*Re+(1-t)*Ri)*cos(a), (t*Re+(1-t)*Ri)*sin(a), z];
            end

            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, fens.xyz,
            #     FinEtools.MeshExportModule.H8)
            # @async run(`"paraview.exe" $File`)

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

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            # Plane-stress constraint: assume the plane z=0 is the plane of symmetry of the plate
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)
            # If this was enabled, the plane-strain  constraint would be enforced.
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)
            el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), GaussRule(2, 2)))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                local r = sqrt(XYZ[1]^2 + XYZ[2]^2)
                nx = XYZ[1]/r; ny = XYZ[2]/r
                # local sx, sy, txy
                # sx, sy, txy = sigmaxx(XYZ), sigmayy(XYZ), sigmaxy(XYZ)
                # sn = sx * nx^2 + sy * ny^2 + 2 * nx * ny * txy
                # tn = -(sx - sy) * nx * ny + (nx^2 - ny^2) * txy
                # forceout[1] = sn * nx - tn * ny
                # forceout[2] = sn * ny + tn * nx
                # forceout[3] = 0.0
                forceout[1] = sigmarr(XYZ) * nx - sigmart(XYZ) * ny
                forceout[2] = sigmarr(XYZ) * ny + sigmart(XYZ) * nx
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

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

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

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

            nlA = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            nlB = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, Thickness], inflate=tolerance);
            # thecorneru = zeros(FFlt,length(nlA),3)
            # gathervalues_asmat!(u, thecorneru, nl);
            # thecorneru = mean(thecorneru, 1)[1]/phun("mm")
            # println("displacement = $(thecorneru) vs -0.10215 [MM]")

            println("Extrapolation: $( extrapolation )")
            sigx = fieldfromintegpoints(femm, geom, u, :Cauchy, 1;
            nodevalmethod = :averaging, reportat = extrapolation)
            sigy = fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
            nodevalmethod = :averaging, reportat = extrapolation)
            sigyA = mean(sigy.values[nlA,1], dims = 1)[1]
            sigyAtrue = sigmatt([Ri, 0.0, 0.0])
            println("sig_y@A =$(sigyA/phun("MPa")) vs $(sigyAtrue/phun("MPa")) [MPa]")
            sigxB = mean(sigx.values[nlB,1], dims = 1)[1]
            sigxBtrue = sigmatt([0.0, Ri, 0.0])
            println("sig_x@B =$(sigxB/phun("MPa")) vs $(sigxBtrue/phun("MPa")) [MPa]")
            push!(numnodes, count(fens))
            push!(numelements, count(fes))
            push!(sigxderrs[extrapolation], abs(sigxB/sigxBtrue - 1.0))
            push!(sigyderrs[extrapolation], abs(sigyA/sigyAtrue - 1.0))
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            # FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
            # scalars=[("sigmax", sigx.values/phun("MEGA*PA"))])
            # @async run(`"paraview.exe" $File`)
        end
    end

    File = "plate_w_hole_PE_MSH8_convergence.CSV"
    savecsv(File, numelements=vec(numelements), numnodes=vec(numnodes), sigxderrtrend=vec(sigxderrs[:extraptrend]), sigxderrdefault=vec(sigxderrs[:extrapmean]), sigyderrtrend=vec(sigyderrs[:extraptrend]), sigyderrdefault=vec(sigyderrs[:extrapmean]))

end # plate_w_hole_MSH8_PE_convergence

function plate_w_hole_MST10_convergence()
    E = 2.4*phun("MEGA*PA");# 210e3 MPa
    nu = 0.49995;
    Re = 0.3*phun("M"); # outer radius
    Ri= 0.1*phun("M"); # hole radius
    H = 0.1*phun("M") # thickness of the plate
    nRadial, nCircumferential=3, 5;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    sigxderrs = Dict{Symbol, FFltVec}()
    sigyderrs = Dict{Symbol, FFltVec}()
    numelements = []
    numnodes = []
    for extrapolation in [:extraptrend :extrapmean]
        sigxderrs[extrapolation] = FFltVec[]
        sigyderrs[extrapolation] = FFltVec[]
        numelements = []
        numnodes = []
        for ref in 0:1:3
            # Thickness = H
            Thickness = H/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = T10block(1.0, pi/2, Thickness, 2^ref*nRadial, 2^ref*nCircumferential, 1)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);

            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
                fens.xyz[i,:] = [(t*Re+(1-t)*Ri)*cos(a), (t*Re+(1-t)*Ri)*sin(a), z];
            end

            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, fens.xyz,
            #     FinEtools.MeshExportModule.H8)
            # @async run(`"paraview.exe" $File`)

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

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            # Plane-stress constraint: assume the plane z=0 is the plane of symmetry of the plate
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)
            # If this was enabled, the plane-strain  constraint would be enforced.
            # l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
            # setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)
            el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), TriRule(3)))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                local r = sqrt(XYZ[1]^2 + XYZ[2]^2)
                nx = XYZ[1]/r; ny = XYZ[2]/r
                # local sx, sy, txy
                # sx, sy, txy = sigmaxx(XYZ), sigmayy(XYZ), sigmaxy(XYZ)
                # sn = sx * nx^2 + sy * ny^2 + 2 * nx * ny * txy
                # tn = -(sx - sy) * nx * ny + (nx^2 - ny^2) * txy
                # forceout[1] = sn * nx - tn * ny
                # forceout[2] = sn * ny + tn * nx
                # forceout[3] = 0.0
                forceout[1] = sigmarr(XYZ) * nx - sigmart(XYZ) * ny
                forceout[2] = sigmarr(XYZ) * ny + sigmart(XYZ) * nx
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinearMST10(MR, IntegDomain(fes, TetRule(4)), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

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

            nlA = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            nlB = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, Thickness], inflate=tolerance);
            # thecorneru = zeros(FFlt,length(nlA),3)
            # gathervalues_asmat!(u, thecorneru, nl);
            # thecorneru = mean(thecorneru, 1)[1]/phun("mm")
            # println("displacement = $(thecorneru) vs -0.10215 [MM]")

            println("Extrapolation: $( extrapolation )")
            sigx = fieldfromintegpoints(femm, geom, u, :Cauchy, 1;            nodevalmethod = :averaging, reportat = extrapolation)
            sigy = fieldfromintegpoints(femm, geom, u, :Cauchy, 2;            nodevalmethod = :averaging, reportat = extrapolation)
            sigyA = mean(sigy.values[nlA,1], dims = 1)[1]
            sigyAtrue = sigmatt([Ri, 0.0, 0.0])
            println("sig_y@A =$(sigyA/phun("MPa")) vs $(sigyAtrue/phun("MPa")) [MPa]")
            sigxB = mean(sigx.values[nlB,1], dims = 1)[1]
            sigxBtrue = sigmatt([0.0, Ri, 0.0])
            println("sig_x@B =$(sigxB/phun("MPa")) vs $(sigxBtrue/phun("MPa")) [MPa]")
            push!(numnodes, count(fens))
            push!(numelements, count(fes))
            push!(sigxderrs[extrapolation], abs(sigxB/sigxBtrue - 1.0))
            push!(sigyderrs[extrapolation], abs(sigyA/sigyAtrue - 1.0))
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            # FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
            # scalars=[("sigmax", sigx.values/phun("MEGA*PA"))])
            # @async run(`"paraview.exe" $File`)
        end
    end

    File = "plate_w_hole_MST10_convergence.CSV"
    savecsv(File, numelements=vec(numelements), numnodes=vec(numnodes), sigxderrtrend=vec(sigxderrs[:extraptrend]), sigxderrdefault=vec(sigxderrs[:extrapmean]), sigyderrtrend=vec(sigyderrs[:extraptrend]), sigyderrdefault=vec(sigyderrs[:extrapmean]))

end # plate_w_hole_MST10_convergence

function plate_w_hole_MST10_PE_convergence()
    E = 2.4*phun("MEGA*PA");# 210e3 MPa
    nu = 0.49995;
    Re = 0.3*phun("M"); # outer radius
    Ri= 0.1*phun("M"); # hole radius
    H = 0.1*phun("M") # thickness of the plate
    nRadial, nCircumferential=3, 5;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    sigxderrs = Dict{Symbol, FFltVec}()
    sigyderrs = Dict{Symbol, FFltVec}()
    numelements = []
    numnodes = []
    for extrapolation in [:extraptrend :extrapmean]
        sigxderrs[extrapolation] = FFltVec[]
        sigyderrs[extrapolation] = FFltVec[]
        numelements = []
        numnodes = []
        for ref in 0:1:3
            # Thickness = H
            Thickness = H/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = T10block(1.0, pi/2, Thickness, 2^ref*nRadial, 2^ref*nCircumferential, 1)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);

            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
                fens.xyz[i,:] = [(t*Re+(1-t)*Ri)*cos(a), (t*Re+(1-t)*Ri)*sin(a), z];
            end

            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, fens.xyz,
            #     FinEtools.MeshExportModule.H8)
            # @async run(`"paraview.exe" $File`)

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

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            # Plane-stress constraint: assume the plane z=0 is the plane of symmetry of the plate
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)
            # If this was enabled, the plane-strain  constraint would be enforced.
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)
            el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), TriRule(3)))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                local r = sqrt(XYZ[1]^2 + XYZ[2]^2)
                nx = XYZ[1]/r; ny = XYZ[2]/r
                # local sx, sy, txy
                # sx, sy, txy = sigmaxx(XYZ), sigmayy(XYZ), sigmaxy(XYZ)
                # sn = sx * nx^2 + sy * ny^2 + 2 * nx * ny * txy
                # tn = -(sx - sy) * nx * ny + (nx^2 - ny^2) * txy
                # forceout[1] = sn * nx - tn * ny
                # forceout[2] = sn * ny + tn * nx
                # forceout[3] = 0.0
                forceout[1] = sigmarr(XYZ) * nx - sigmart(XYZ) * ny
                forceout[2] = sigmarr(XYZ) * ny + sigmart(XYZ) * nx
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinearMST10(MR, IntegDomain(fes, TetRule(4)), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

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

            nlA = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            nlB = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, Thickness], inflate=tolerance);
            # thecorneru = zeros(FFlt,length(nlA),3)
            # gathervalues_asmat!(u, thecorneru, nl);
            # thecorneru = mean(thecorneru, 1)[1]/phun("mm")
            # println("displacement = $(thecorneru) vs -0.10215 [MM]")

            println("Extrapolation: $( extrapolation )")
            sigx = fieldfromintegpoints(femm, geom, u, :Cauchy, 1;
            nodevalmethod = :averaging, reportat = extrapolation)
            sigy = fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
            nodevalmethod = :averaging, reportat = extrapolation)
            sigyA = mean(sigy.values[nlA,1], dims = 1)[1]
            sigyAtrue = sigmatt([Ri, 0.0, 0.0])
            println("sig_y@A =$(sigyA/phun("MPa")) vs $(sigyAtrue/phun("MPa")) [MPa]")
            sigxB = mean(sigx.values[nlB,1], dims = 1)[1]
            sigxBtrue = sigmatt([0.0, Ri, 0.0])
            println("sig_x@B =$(sigxB/phun("MPa")) vs $(sigxBtrue/phun("MPa")) [MPa]")
            push!(numnodes, count(fens))
            push!(numelements, count(fes))
            push!(sigxderrs[extrapolation], abs(sigxB/sigxBtrue - 1.0))
            push!(sigyderrs[extrapolation], abs(sigyA/sigyAtrue - 1.0))
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            # FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
            # scalars=[("sigmax", sigx.values/phun("MEGA*PA"))])
            # @async run(`"paraview.exe" $File`)
        end
    end

    File = "plate_w_hole_PE_MST10_convergence.CSV"
    savecsv(File, numelements=vec(numelements), numnodes=vec(numnodes), sigxderrtrend=vec(sigxderrs[:extraptrend]), sigxderrdefault=vec(sigxderrs[:extrapmean]), sigyderrtrend=vec(sigyderrs[:extraptrend]), sigyderrdefault=vec(sigyderrs[:extrapmean]))

end # plate_w_hole_MST10_PE_convergence

function plate_w_hole_MST10_stress()
    E = 2.4*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    Re = 0.3*phun("M"); # outer radius
    Ri= 0.1*phun("M"); # hole radius
    H = 0.1*phun("M") # thickness of the plate
    nRadial, nCircumferential, nThickness = 6, 8, 1;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    convergencestudy = FDataDict[]
    for ref in 0:1:2
        println("ref = $(ref)")
        # Thickness = H
        Thickness = H/2^ref
        tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

        fens,fes = T10block(1.0, pi/2, Thickness, 2^ref*nRadial, 2^ref*nCircumferential, 2^ref*nThickness)

        bdryfes = meshboundary(fes);
        icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);

        for i=1:count(fens)
            t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
            fens.xyz[i,:] = [(t*Re+(1-t)*Ri)*cos(a), (t*Re+(1-t)*Ri)*sin(a), z];
        end

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

        l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
        setebc!(u,l1,true, 2, 0.0)
        l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
        setebc!(u,l1,true, 1, 0.0)
        # Plane-stress constraint: assume the plane z=0 is the plane of symmetry of the plate
        l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
        setebc!(u,l1,true, 3, 0.0)
        # If this was enabled, the plane-strain  constraint would be enforced.
        # l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
        # setebc!(u,l1,true, 3, 0.0)

        applyebc!(u)
        numberdofs!(u)
        el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), TriRule(3)))
        function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
            local r = sqrt(XYZ[1]^2 + XYZ[2]^2)
            nx = XYZ[1]/r; ny = XYZ[2]/r
            forceout[1] = sigmarr(XYZ) * nx - sigmart(XYZ) * ny
            forceout[2] = sigmarr(XYZ) * ny + sigmart(XYZ) * nx
            forceout[3] = 0.0
            return forceout
        end
        fi = ForceIntensity(FFlt, 3, pfun);
        F2 = distribloads(el1femm, geom, u, fi, 2);

        MR = DeforModelRed3D

        material = MatDeforElastIso(MR, E, nu)

        femm = FEMMDeforLinearMST10(MR, IntegDomain(fes, TetRule(4)), material)

        # The geometry field now needs to be associated with the FEMM
        femm = associategeometry!(femm, geom)

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

        stressfields = elemfieldfromintegpoints(femm, geom, u, :Cauchy, collect(1:6))

        push!(convergencestudy, FDataDict(
        "elementsize"=> 1.0 / 2^ref,
        "fens"=>fens,
        "fes"=>fes,
        "geom"=>geom,
        "u"=>u,
        "femm"=>femm,
        "stressfields"=>[stressfields],
        "tolerance"=>tolerance)
        )
    end # for ref in

    # File = "mplate_w_hole_MST10m_stress"
    # open(File * ".jls", "w") do file
    #     serialize(file, convergencestudy)
    # end

end # plate_w_hole_MST10_stress

function plate_w_hole_RECT_H20_convergence()
    E = 210000*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    Ri= 0.15*phun("M"); # hole radius
    Re = 2*Ri; # outer radius
    H = 0.01*phun("M") # thickness of the plate
    nRadial, nCircumferential=6, 3;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(vec(x[1:2]));
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    sigyderrs = Dict{Symbol, FFltVec}()

    nelems = []
    for extrapolation in [:meanonly]
        sigyderrs[extrapolation] = FFltVec[]
        nelems = []
        for ref in [1]
            Thickness = H
            # Thickness = H/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = H8elliphole(Ri, Ri, Re, Re, Thickness,
            2^ref*nCircumferential, 2^ref*nCircumferential, 2^ref*nRadial, 1)
            fens,fes = H8toH20(fens,fes)
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, fens.xyz,
            #     FinEtools.MeshExportModule.H20)
            # @async run(`"paraview.cexe" $File`)

            println("My mesh=>$((count(fens), count(fes)))")
            #
            # output = import_ABAQUS("plane_w_hole_m_debug.inp")
            # fens1,fes1 = output["fens"], output["fesets"][1]
            # println("Matlab mesh=>$((count(fens1), count(fes1[1])))")
            #
            #  fens3, newfes1, fes2 = mergemeshes(fens,fes, fens1,fes1[1], tolerance)
            #  fes3 = cat(2, newfes1)
            #  println("Merged mesh=>$((count(fens3), count(fes3)))")

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

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)
            # l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
            # setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)

            bdryfes = meshboundary(fes);
            # ixl = selectelem(fens, bdryfes, plane=[1.0, 0.0, 0.0, Re], thickness=tolerance);
            ixl = selectelem(fens, bdryfes, box=[Re, Re, -Inf, +Inf, -Inf, +Inf], inflate = tolerance);
            elxfemm =  FEMMBase(IntegDomain(subset(bdryfes,ixl), GaussRule(2, 2)))
            function pfunx(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                forceout[1] = sigmaxx(XYZ)
                forceout[2] = sigmaxy(XYZ)
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfunx);
            Fx = distribloads(elxfemm, geom, u, fi, 2);
            # iyl = selectelem(fens, bdryfes, plane=[0.0, 1.0, 0.0, Re], thickness=tolerance);
            iyl = selectelem(fens, bdryfes, box=[-Inf, +Inf, Re, Re, -Inf, +Inf], inflate = tolerance);
            elyfemm =  FEMMBase(IntegDomain(subset(bdryfes,iyl), GaussRule(2, 2)))
            function pfuny(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                forceout[1] = sigmaxy(XYZ)
                forceout[2] = sigmayy(XYZ)
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfuny);
            Fy = distribloads(elyfemm, geom, u, fi, 2);

            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

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

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

            K = stiffness(femm, geom, u)
            K = cholesky(K)
            U = K\(Fx + Fy)
            scattersysvec!(u, U[:])
            println("oof load = $(norm(Fx + Fy, 2))")

            nlA = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, 00.0], inflate=tolerance);
            pointu = zeros(FFlt,length(nlA),3)
            gathervalues_asmat!(u, pointu, nlA);
            println("disp@A = $(pointu/phun("mm")) [MM]")
            nlB = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, 0.0], inflate=tolerance);
            pointu = zeros(FFlt,length(nlB),3)
            gathervalues_asmat!(u, pointu, nlB);
            println("disp@B = $(pointu/phun("mm")) [MM]")
            nlC = selectnode(fens, box=[Re, Re, Re, Re, Thickness, Thickness], inflate=tolerance);
            pointu = zeros(FFlt,length(nlC),3)
            gathervalues_asmat!(u, pointu, nlC);
            println("disp@C = $(pointu/phun("mm")) [MM]")

            nlAallz = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            nlBallz = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, Thickness], inflate=tolerance);
            sigx = fieldfromintegpoints(femm, geom, u, :Cauchy, 1; nodevalmethod = :invdistance, reportat = extrapolation)
            sigy = fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :invdistance, reportat = extrapolation)
            sigyA = mean(sigy.values[nlAallz,1], dims = 1)[1]
            sigyAtrue = sigmayy([Ri, 0.0, 0.0])
            println("sig_y@A =$(sigyA/phun("MPa")) vs $(sigyAtrue/phun("MPa")) [MPa]")
            sigxB = mean(sigx.values[nlBallz,1], dims = 1)[1]
            sigxBtrue = sigmaxx([0.0, Ri, 0.0])
            println("sig_x@B =$(sigxB/phun("MPa")) vs $(sigxBtrue/phun("MPa")) [MPa]")
            # println("$extrapolation, $(count(fes)), $(sigyd/phun("MPa"))")
            # push!(nelems, count(fes))
            # push!(sigyderrs[extrapolation], abs(sigyd/sigma_yD - 1.0))
            File =  "a.vtk"
            vtkexportmesh(File, fes.conn, geom.values,
            FinEtools.MeshExportModule.H20; vectors=[("u", u.values)],
            scalars=[("sigmax", sigx.values/phun("MEGA*PA")),
            ("sigmay", sigy.values/phun("MEGA*PA"))])
            @async run(`"paraview.exe" $File`)
        end
    end

    # df = DataFrame(nelems=vec(nelems),
    #     sigyderrtrend=vec(sigyderrs[:extraptrend]),
    #     sigyderrdefault=vec(sigyderrs[:extrapmean]))
    # File = "LE1NAFEMS_MSH8_convergence.CSV"
    # CSV.write(File, df)
    # @async run(`"paraview.exe" $File`)

end # plate_w_hole_RECT_H20_convergence

function plate_w_hole_RECT_MSH8_convergence()
    E = 210000*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    Ri= 0.1*phun("M"); # hole radius
    Re = 2*Ri; # outer radius
    H = 0.01*phun("M") # thickness of the plate
    nRadial, nCircumferential=6, 3;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    sigyderrs = Dict{Symbol, FFltVec}()

    nelems = []
    for extrapolation in [:extraptrend :extrapmean]
        sigyderrs[extrapolation] = FFltVec[]
        nelems = []
        for ref in [2]
            # Thickness = H
            Thickness = H/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = H8elliphole(Ri, Ri, Re, Re, Thickness,
            2^ref*nCircumferential, 2^ref*nCircumferential, 2^ref*nRadial, 1)

            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, fens.xyz,
            #     FinEtools.MeshExportModule.H8)
            # @async run(`"paraview.exe" $File`)

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

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)
            # l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
            # setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)

            bdryfes = meshboundary(fes);
            ixl = selectelem(fens, bdryfes, plane=[1.0, 0.0, 0.0, Re], thickness=tolerance);
            elxfemm =  FEMMBase(IntegDomain(subset(bdryfes,ixl), GaussRule(2, 2)))
            function pfunx(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                forceout[1] = sigmaxx(XYZ)
                forceout[2] = sigmaxy(XYZ)
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfunx);
            Fx = distribloads(elxfemm, geom, u, fi, 2);
            iyl = selectelem(fens, bdryfes, plane=[0.0, 1.0, 0.0, Re], thickness=tolerance);
            elyfemm =  FEMMBase(IntegDomain(subset(bdryfes,iyl), GaussRule(2, 2)))
            function pfuny(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                forceout[1] = -sigmaxy(XYZ)
                forceout[2] = sigmayy(XYZ)
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfuny);
            Fy = distribloads(elyfemm, geom, u, fi, 2);

            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

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

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

            K = stiffness(femm, geom, u)
            K = cholesky(K)
            U = K\(Fx + Fy)
            scattersysvec!(u, U[:])

            nlA = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            nlB = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, Thickness], inflate=tolerance);
            # thecorneru = zeros(FFlt,length(nlA),3)
            # gathervalues_asmat!(u, thecorneru, nl);
            # thecorneru = mean(thecorneru, 1)[1]/phun("mm")
            # println("displacement = $(thecorneru) vs -0.10215 [MM]")

            println("Extrapolation: $( extrapolation )")
            sigx = fieldfromintegpoints(femm, geom, u, :Cauchy, 1; nodevalmethod = :averaging, reportat = extrapolation)
            sigy = fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :averaging, reportat = extrapolation)
            sigyA = mean(sigy.values[nlA,1], dims = 1)[1]
            sigyAtrue = sigmatt([Ri, 0.0, 0.0])
            println("sig_y@A =$(sigyA/phun("MPa")) vs $(sigyAtrue/phun("MPa")) [MPa]")
            sigxB = mean(sigx.values[nlB,1], dims = 1)[1]
            sigxBtrue = sigmatt([0.0, Ri, 0.0])
            println("sig_x@B =$(sigxB/phun("MPa")) vs $(sigxBtrue/phun("MPa")) [MPa]")
            # println("$extrapolation, $(count(fes)), $(sigyd/phun("MPa"))")
            # push!(nelems, count(fes))
            # push!(sigyderrs[extrapolation], abs(sigyd/sigma_yD - 1.0))
            File =  "a.vtk"
            vtkexportmesh(File, fes.conn, geom.values,
            FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
            scalars=[("sigmax", sigx.values/phun("MEGA*PA")),
            ("sigmay", sigy.values/phun("MEGA*PA"))])
            @async run(`"paraview.exe" $File`)
        end
    end

    # df = DataFrame(nelems=vec(nelems),
    #     sigyderrtrend=vec(sigyderrs[:extraptrend]),
    #     sigyderrdefault=vec(sigyderrs[:extrapmean]))
    # File = "LE1NAFEMS_MSH8_convergence.CSV"
    # CSV.write(File, df)
    # @async run(`"paraview.exe" $File`)

end # plate_w_hole_RECT_MSH8_convergence

function plate_w_hole_T10_PE_convergence()
    E = 2.4*phun("MEGA*PA");# 210e3 MPa
    nu = 0.49995;
    Re = 0.3*phun("M"); # outer radius
    Ri= 0.1*phun("M"); # hole radius
    H = 0.1*phun("M") # thickness of the plate
    nRadial, nCircumferential=3, 5;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    sigxderrs = Dict{Symbol, FFltVec}()
    sigyderrs = Dict{Symbol, FFltVec}()
    numelements = []
    numnodes = []
    for extrapolation in [:extraptrend :extrapmean]
        sigxderrs[extrapolation] = FFltVec[]
        sigyderrs[extrapolation] = FFltVec[]
        numelements = []
        numnodes = []
        for ref in 0:1:2
            # Thickness = H
            Thickness = H/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = T10block(1.0, pi/2, Thickness, 2^ref*nRadial, 2^ref*nCircumferential, 1)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);

            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
                fens.xyz[i,:] = [(t*Re+(1-t)*Ri)*cos(a), (t*Re+(1-t)*Ri)*sin(a), z];
            end

            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, fens.xyz,
            #     FinEtools.MeshExportModule.H8)
            # @async run(`"paraview.exe" $File`)

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

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            # Plane-stress constraint: assume the plane z=0 is the plane of symmetry of the plate
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)
            # If this was enabled, the plane-strain  constraint would be enforced.
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)
            el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), TriRule(3)))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                local r = sqrt(XYZ[1]^2 + XYZ[2]^2)
                nx = XYZ[1]/r; ny = XYZ[2]/r
                # local sx, sy, txy
                # sx, sy, txy = sigmaxx(XYZ), sigmayy(XYZ), sigmaxy(XYZ)
                # sn = sx * nx^2 + sy * ny^2 + 2 * nx * ny * txy
                # tn = -(sx - sy) * nx * ny + (nx^2 - ny^2) * txy
                # forceout[1] = sn * nx - tn * ny
                # forceout[2] = sn * ny + tn * nx
                # forceout[3] = 0.0
                forceout[1] = sigmarr(XYZ) * nx - sigmart(XYZ) * ny
                forceout[2] = sigmarr(XYZ) * ny + sigmart(XYZ) * nx
                forceout[3] = 0.0
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinear(MR, IntegDomain(fes, TetRule(4)), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

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

            nlA = selectnode(fens, box=[Ri, Ri, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            nlB = selectnode(fens, box=[0.0, 0.0, Ri, Ri, 0.0, Thickness], inflate=tolerance);
            # thecorneru = zeros(FFlt,length(nlA),3)
            # gathervalues_asmat!(u, thecorneru, nl);
            # thecorneru = mean(thecorneru, 1)[1]/phun("mm")
            # println("displacement = $(thecorneru) vs -0.10215 [MM]")

            println("Extrapolation: $( extrapolation )")
            sigx = fieldfromintegpoints(femm, geom, u, :Cauchy, 1)
            sigy = fieldfromintegpoints(femm, geom, u, :Cauchy, 2)
            sigyA = mean(sigy.values[nlA,1], dims = 1)[1]
            sigyAtrue = sigmatt([Ri, 0.0, 0.0])
            println("sig_y@A =$(sigyA/phun("MPa")) vs $(sigyAtrue/phun("MPa")) [MPa]")
            sigxB = mean(sigx.values[nlB,1], dims = 1)[1]
            sigxBtrue = sigmatt([0.0, Ri, 0.0])
            println("sig_x@B =$(sigxB/phun("MPa")) vs $(sigxBtrue/phun("MPa")) [MPa]")
            push!(numnodes, count(fens))
            push!(numelements, count(fes))
            push!(sigxderrs[extrapolation], abs(sigxB/sigxBtrue - 1.0))
            push!(sigyderrs[extrapolation], abs(sigyA/sigyAtrue - 1.0))
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            # FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
            # scalars=[("sigmax", sigx.values/phun("MEGA*PA"))])
            # @async run(`"paraview.exe" $File`)
        end
    end

    File = "plate_w_hole_PE_T10_convergence.CSV"
    savecsv(File, numelements=vec(numelements), numnodes=vec(numnodes), sigxderrtrend=vec(sigxderrs[:extraptrend]), sigxderrdefault=vec(sigxderrs[:extrapmean]), sigyderrtrend=vec(sigyderrs[:extraptrend]), sigyderrdefault=vec(sigyderrs[:extrapmean]))

end # plate_w_hole_T10_PE_convergence

function plate_w_hole_T10_stress()
    E = 2.4*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    Re = 0.3*phun("M"); # outer radius
    Ri= 0.1*phun("M"); # hole radius
    H = 0.1*phun("M") # thickness of the plate
    nRadial, nCircumferential, nThickness = 6, 8, 1;
    sigma0=1*phun("MEGA*PA");

    function sigmaxx(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0*(1-Ri^2/r^2*(3/2*cos(2*th)+cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmayy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*cos(2*th)-cos(4*th))+3/2*Ri^4/r^4*cos(4*th));
    end
    function sigmaxy(x)
        local r=norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0*(Ri^2/r^2*(1/2*sin(2*th)+sin(4*th))-3/2*Ri^4/r^4*sin(4*th));
    end
    function sigmarr(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1-Ri^2/r^2) + sigma0/2*(1-4*Ri^2/r^2+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmatt(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return sigma0/2*(1+Ri^2/r^2) - sigma0/2*(1+3*Ri^4/r^4)*cos(2*th)
    end
    function sigmart(x)
        local r = norm(x[1:2]);
        local th = atan(x[2],x[1]);
        return -sigma0/2*(1+2*Ri^2/r^2-3*Ri^4/r^4)*sin(2*th)
    end

    convergencestudy = FDataDict[]
    for ref in 0:1:2
        println("ref = $(ref)")
        # Thickness = H
        Thickness = H/2^ref
        tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

        fens,fes = T10block(1.0, pi/2, Thickness, 2^ref*nRadial, 2^ref*nCircumferential, 2^ref*nThickness)

        bdryfes = meshboundary(fes);
        icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);

        for i=1:count(fens)
            t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
            fens.xyz[i,:] = [(t*Re+(1-t)*Ri)*cos(a), (t*Re+(1-t)*Ri)*sin(a), z];
        end

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

        l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
        setebc!(u,l1,true, 2, 0.0)
        l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
        setebc!(u,l1,true, 1, 0.0)
        # Plane-stress constraint: assume the plane z=0 is the plane of symmetry of the plate
        l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
        setebc!(u,l1,true, 3, 0.0)
        # If this was enabled, the plane-strain  constraint would be enforced.
        # l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, Thickness, Thickness], inflate = tolerance)
        # setebc!(u,l1,true, 3, 0.0)

        applyebc!(u)
        numberdofs!(u)
        el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), TriRule(3)))
        function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
            local r = sqrt(XYZ[1]^2 + XYZ[2]^2)
            nx = XYZ[1]/r; ny = XYZ[2]/r
            forceout[1] = sigmarr(XYZ) * nx - sigmart(XYZ) * ny
            forceout[2] = sigmarr(XYZ) * ny + sigmart(XYZ) * nx
            forceout[3] = 0.0
            return forceout
        end
        fi = ForceIntensity(FFlt, 3, pfun);
        F2 = distribloads(el1femm, geom, u, fi, 2);

        MR = DeforModelRed3D

        material = MatDeforElastIso(MR, E, nu)

        femm = FEMMDeforLinear(MR, IntegDomain(fes, TetRule(4)), material)

        # The geometry field now needs to be associated with the FEMM
        femm = associategeometry!(femm, geom)

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

        stressfields = elemfieldfromintegpoints(femm, geom, u, :Cauchy, collect(1:6))

        push!(convergencestudy, FDataDict(
        "elementsize"=> 1.0 / 2^ref,
        "fens"=>fens,
        "fes"=>fes,
        "geom"=>geom,
        "u"=>u,
        "femm"=>femm,
        "stressfields"=>[stressfields],
        "tolerance"=>tolerance)
        )
    end # for ref in

    # File = "mplate_w_hole_T10m_stress"
    # open(File * ".jls", "w") do file
    #     serialize(file, convergencestudy)
    # end

end # plate_w_hole_T10_stress

function allrun()
    println("#####################################################")
    println("# plate_w_hole_T10_PE_convergence ")
    plate_w_hole_T10_PE_convergence()
    # println("#####################################################")
    # println("# plate_w_hole_T10_stress ")
    # plate_w_hole_T10_stress()
    # println("#####################################################")
    # println("# plate_w_hole_RECT_MSH8_convergence ")
    # plate_w_hole_RECT_MSH8_convergence()
    # println("#####################################################")
    # println("# plate_w_hole_H20_stress ")
    # plate_w_hole_H20_stress()
    # println("#####################################################")
    # println("# plate_w_hole_MSH8_convergence ")
    # plate_w_hole_MSH8_convergence()
    # println("#####################################################")
    # println("# plate_w_hole_MSH8_PE_convergence ")
    # plate_w_hole_MSH8_PE_convergence()
    # println("#####################################################")
    # println("# plate_w_hole_MST10_convergence ")
    # plate_w_hole_MST10_convergence()
    # println("#####################################################")
    # println("# plate_w_hole_MST10_PE_convergence ")
    # plate_w_hole_MST10_PE_convergence()
    # println("#####################################################")
    # println("# plate_w_hole_MST10_stress ")
    # plate_w_hole_MST10_stress()
    # println("#####################################################")
    # println("# plate_w_hole_RECT_H20_convergence ")
    # plate_w_hole_RECT_H20_convergence()
    return true
end # function allrun

end # module plate_w_hole_examples


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Model of basic streams
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt and Rafael de P. Soares
#* $Id$
#*---------------------------------------------------------------------
type streamPHS
	streamPHS()=begin
		new(
			streamPH(),
			entr_mol(Dict{Symbol,Any}(
				:Brief=>"Stream Entropy",
				:Protected=>true
			)),
			[
				:(s = (1-v)*_base_1.PP.LiquidEntropy(T, P, x) + v*_base_1.PP.VapourEntropy(T, P, y)),
			],
			[
				"Entropy",
			],
			[:s,]
		)
	end
	_base_1::streamPH
	s::entr_mol
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export streamPHS
function setEquationFlow(in::streamPHS)
	addEquation(1)
end
function atributes(in::streamPHS,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Brief]="Stream with built-in flash calculation"
	fields[:Info]="
	This model should be used when the vaporization fraction
	is unknown.

	The built-in flash calculation will determine the stream
	state as a function of the overall composition '''z''', the
	pressure '''P''' and the enthalpy '''h'''.

	Additionally, the liquid composition '''x''', the vapor
	composition '''y''' and the stream entropy are calculated.
	"
	fields[:Pallete]=false
	drive!(fields,_)
	return fields
end
streamPHS(_::Dict{Symbol,Any})=begin
	newModel=streamPHS()
	newModel.attributes=atributes(newModel,_)
	newModel
end

ext: jl


text:
# include("simulate.jl")
using FiniteDifferences
using Ipopt
using GenericLinearAlgebra
fdm = central_fdm(5,1);

## Ipopt directly
function trim_FinDiff(j::Vector{Joint},GravityInInertial::Vector{Float64};
  ix::Vector{Integer}=zeros(Integer,20*(nB-1)+1),
  iy::Vector{Integer}=zeros(Integer,7*(nB-1)))
  nB = length(j) + 1; # number of bodies in the joint tree
  # ix: 0 indicates that the corresponding variable is free during optimization

  # iy: 0 indicates that the corresponding second derivative term needs to be zero
  # (specify indices for iy starting from second body)
  # iy = [1;4;5] means that the ẍ. β̈1 and β̈2 are free for j[1].RB2

  # Find which variables are free for optimization
  freeVarInd = findall(x->x==0,ix); # indices of free variables

  # find which acceleration-level terms are desired to be 0/constrained
  freeAccInd = findall(x->x==1,iy) # indices of free accelerations

  # Initial guesses
  # Initial guess is provided initial location
  x0 = []; # Initial guess
  for k = 1:length(j)
    append!(x0,j[k].RB2.x)
  end
  extFListCopy = externalFTotal(0.0,j);
  u0 = genU(extFListCopy); # Initial guess for control input

  n_xu = length(freeVarInd)-1; # Number of optimization variables
  lb = [-Inf*ones(n_xu);0.0]; ub = [Inf*ones(n_xu);Inf] # Bounds on optimization variables
  # lb[4] = 0.0;
  # lb[14+4] = 0.0;

  m = 1; # Number of nonlinear constraints (xdot = f(x,u)) for second order variables
  lb_g = [-Inf]; ub_g = [0.0];

  # extract states and control from optimization variables
  function getXU(z::AbstractArray{T}) where T<:Real
    z0 = [x0;u0[:,2:end][:];0.0];
    z0[freeVarInd] .= 0;

    Mz = zeros(length(z0), length(freeVarInd))
    for (i,ind) in enumerate(freeVarInd)
      Mz[ind,i] = 1;
    end

    zFull = z0 + Mz*z;

    X = zFull[1:14*(nB-1)];
    U = zFull[14*(nB-1)+1:end-1]
    return X,U
  end

  ## Constraint functions
  function qddFn(xuVar)
    x = xuVar[1:14*nB]; u = reshape(xuVar[14*nB+1:end],(6,nB));
    xdot = fxdot(x,u,j,g);
    y = Vector(undef,7*nB); # qdd
    for i=1:nB
      y[7*(i-1)+1:7*(i)] = xdot[7*(2*i-1)+1:14*i]
    end
    deleteat!(y,1:7); # inertial frame has no acceleration
    deleteat!(y,freeAccInd); # free acceleration variables
    return y
  end

  function quatNormConstr(xuVar)
    x = xuVar[1:14*nB]; u = reshape(xuVar[14*nB+1:end],(6,nB));
    q_constr = Vector(undef, nB)
    for i=1:nB
        q_constr[i] = norm(x[14*(i-1)+4:14*(i-1)+7])^2 - 1;
    end
    return q_constr
  end

  function jointLocConstr(xuVar)
    x = xuVar[1:14*nB]; u = reshape(xuVar[14*nB+1:end],(6,nB));
    jLoc_constr = Matrix(undef,3,length(j))
    for i=1:length(j)
        b1_id = j[i].RB1.bodyID; b2_id = j[i].RB2.bodyID;
        x1 = x[(b1_id-1)*14+1:b1_id*14]
        x2 = x[(b2_id-1)*14+1:b2_id*14]
        j_b1 = x1[1:3] + transpose(quat2dcm(x1[4:7]))*j[i].pos1;
        j_b2 = x2[1:3] + transpose(quat2dcm(x2[4:7]))*j[i].pos2;
        jLoc_constr[:,i] = j_b1 - j_b2;
    end
    return jLoc_constr
  end

  function _nonlinear_constraints(z::AbstractArray{T}) where T<:Real
    xVar, uVar = getXU(z)
    x = [Vector{T}(j[1].RB1.x); xVar]
    u = [zeros(T,6) reshape(uVar,(6, nB-1))]
    xuVar = [x;u[:]]; # Vector of x,u

    # Dynamics constraint
    y = qddFn(xuVar);
    # Quaternion norm constraint
    q_constr = quatNormConstr(xuVar);
    # JointLocation constraint  (translation constraint)
    jLoc_constr = jointLocConstr(xuVar);

    constr = [y;q_constr;jLoc_constr[:]];

    J = sum(constr.*constr) - z[end];
    return J
  end

  function eval_f(z)
    J = z[end]
    return J
  end

  function eval_g(z, g)
    g[:] .= _nonlinear_constraints(z)
  end

  function eval_grad_f(z, grad_f)
    grad_f[:] = FiniteDifferences.grad(fdm, eval_f, z)[1];
  end

  function eval_jac_g(z, mode, rows, cols, values)
    nvar = length(z);
    ng = 1;

    if mode == :Structure
        index = 1;
        for i=1:ng
            for j=1:nvar
                rows[index] = i;
                cols[index] = j;
                index += 1;
            end
        end
    else
        values[:] = FiniteDifferences.grad(fdm,_nonlinear_constraints,z)[1];
    end

    # return jac_g
  end

  function grad_f_fn(z)
    grad_f = FiniteDifferences.grad(fdm, eval_f, z)[1];
    return grad_f;
  end

  function grad_nonlinear_constraints_fn(z)
    grad_nc = FiniteDifferences.grad(fdm, _nonlinear_constraints, z)[1];
    return grad_nc;
  end

  function eval_h(z, mode, rows, cols, obj_factor, lambda, values) # Hessian of cost.
    nvar = length(z);
    if mode == :Structure
      # Symmetric matrix, fill the lower left triangle only
      idx = 1
      for row = 1:nvar
        for col = 1:row
          rows[idx] = row
          cols[idx] = col
          idx += 1
        end
      end
    else
      H = obj_factor*FiniteDifferences.jacobian(fdm, grad_f_fn, z)[1];

      for i=1:length(lambda)
        H += lambda[i]*FiniteDifferences.jacobian(fdm, grad_nonlinear_constraints_fn, z)[1];
      end

      # Copy them into return variable -- only lower triangle.
      idx = 1;
      for row = 1:nvar
        for col = 1:row
          values[idx] = H[row,col];
          idx += 1;
        end
      end
    end
  end

  # Initial guess from freeVarInd
  n = length(freeVarInd);
  guessFull = [x0;u0[:,2:end][:];0.0];
  inGuess = guessFull[freeVarInd];

  nele_jac = Integer(m*n);
  nele_hess = Integer(n*(n+1)/2);
  prob = createProblem(n, lb, ub, m, lb_g, ub_g, nele_jac, nele_hess, eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h);
  prob.x = inGuess; # Initial guess.
  status = solveProblem(prob);

  trim_x = prob.x[1:length(x0)];
  gam = prob.x[end];
  trim_u = prob.x[length(x0)+1:end-1];
  return trim_x, trim_u, gam
end
## testing the trim function
# x0Orig, u0Orig = getXU_0(j);
# nB = length(j) + 1;
# ix= zeros(Integer,20*(nB-1)+1);
# # WITH U=0
# ix[14*(nB-1)+1:20*(nB-1)] .= 1;
# trim_U = u0Orig;
# iy = zeros(Integer,7*(nB-1));
#
# out = trim_FinDiff(j,g,ix=ix, iy=iy);
# trim_x, trim_u, gam = out;
# println("trim_x = ", trim_x)
# trim_X = [j[1].RB1.x;trim_x]
# # trim_U = [zeros(6) reshape(trim_u,(6,nB-1))];
# # println("trim_U = ", trim_U[:,2])
# # trim_U = zeros(6,nB);
# println("qconstr = ", norm(trim_x[4:7])-1)
# println("ẍ =", norm(fxdot(trim_X,trim_U, j, g)))

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
__precompile__()

module PersonalBanking

export main

using ExcelReaders, DataFrames, CSV, StringDistances

include(joinpath(Pkg.dir("PersonalBanking"), "src", "categorise.jl"))

# const rawdata = "/home/yakir/documents/money/rawdata"
const rawdata = "rawdata"

function cleandescription(X)
    x = lowercase(X)
    x = replace(x, r"ö", 'o')
    x = replace(x, r"ä", 'a')
    x = replace(x, r"å", 'a')
    x = replace(x, r"[^a-z0-9]", ' ')
    join(unique(split(x)), ' ')
end

function fixdescription(x)::String
    m = match(r"^[0-9][0-9][0-9][0-9]-[0-9][0-9]-[0-9][0-9]\s?(.*)", x)
    if m ≠ nothing
        strip(first(m.captures))
    else
        strip(x)
    end
end

function excelskandia(file)
    f = openxl(file)
    sheetname = "Kontoutdrag"
    sheet = f.workbook[:sheet_by_name](sheetname)
    startrow, startcol, endrow, endcol = ExcelReaders.convert_args_to_row_col(sheet)
    a = readxl(f, "$(sheetname)!A5:C$endrow")
    description = cleandescription.(fixdescription.(string.(a[:,2])))
    df = DataFrame(date = Date.(strip.(a[:,1])), description = description, amount = float.(a[:,3]))
    sort!(df, :date)
end

fixdate(x) = Date(strip(x)) - Date(0,1,1) + Date(2000,1,1)

function excelswedbank(file)
    f = openxl(file)
    sheetname = "Blad1"
    sheet = f.workbook[:sheet_by_name](sheetname)
    startrow, startcol, endrow, endcol = ExcelReaders.convert_args_to_row_col(sheet)
    date = readxl(f, "$(sheetname)!E6:E$endrow")
    description = readxl(f, "$(sheetname)!G6:G$endrow")
    amount = readxl(f, "$(sheetname)!I6:I$endrow")
    description = cleandescription.(string.(vec(description)))
    df = DataFrame(date = fixdate.(vec(date)), description = description, amount = float.(vec(amount)))
    sort!(df, :date)
end

fixsum(x) = parse(Float64, replace(filter(!isspace, x), ',', '.'))

function csvswedbank(file)
    a = CSV.read(file, delim = ';', header = ["description", "transaktionsdatum", "date", "amount", "saldo"], missingstring = "-")
    dropmissing!(a)
    disallowmissing!(a)
    a[:date] = Date.(a[:date])
    a[:amount] = fixsum.(a[:amount])
    a[:description] .= cleandescription.(a[:description])
    sort!(a[:, [:date, :description, :amount]], :date)
    # a = CSV.read(file, delim = ';', transforms = Dict{Int, Function}(1 => x -> strip(string(x)), 2 => x -> Date(x), 3 => identity, 4 => fixsum, 5 => identity), header = ["description", "a", "date", "amount", "b"])
end

function loadbank(person, ext, file)
    if person == "yakir"
        excelskandia(file)
    else
        if ext == ".csv"
            csvswedbank(file)
        else
            excelswedbank(file)
        end
    end
end

function overlap(t1, t2)
    date, description, amount = convert(Array, t1[end, :])
    date < t2[1, :date] && return 1
    i1 = findfirst(x -> x ≥ date - Dates.Day(1), t2[:date])
    i2 = findlast(x -> x ≤ date + Dates.Day(1), t2[:date])
    i1 == i2 && return i1 + 1
    candidates = filter(i -> amount - 1 ≤ t2[i, :amount] ≤ amount + 1, i1:i2)
    isempty(candidates) && return 1
    length(candidates) == 1 && return candidates[1] + 1
    _, i = findmax(compare(Hamming(), description, t2[candidate, :description]) for candidate in candidates)
    candidates[i] + 1
end

function concatenate!(t1, t2)
    t1, t2 = t1[end, :date] < t2[end, :date] ? (t1, t2) : (t2, t1)
    i = overlap(t1, t2)
    append!(t1, t2[i:end, :])
end

function getall(folder)
    all = Dict("yakir" => Dict{String, DataFrame}(), "ninna" => Dict{String, DataFrame}())
    for f in readdir(folder)
        if f[1] ≠ '.'
            file = joinpath(folder, f)
            person, account, rest = split(f, '_')
            _, ext = splitext(rest)
            bank = loadbank(person, ext, file)
            bank[:account] = Symbol(account)
            bank[:person] = Symbol(person)
            if haskey(all[person], account)
                concatenate!(all[person][account], bank)
            else
                all[person][account] = bank
            end
        end
    end
    all
end

function gather(a)
    table = DataFrame(date = Date(1), description = "", amount = 0.0, account = :a, person = :b)
    deleterows!(table, 1)
    for person in values(a), account in values(person)
        append!(table, account)
    end
    sort!(table, :date)
end

function main(folder)
    a = getall(joinpath(folder, rawdata))
    table = gather(a)
    categorise!(table, folder)
end

end # module

#=badwords = ["malmo", "lund", "till", "montreal", "koebenhavn", "haifa", "fran", "milano", "francisco", "san", "kastrup", "stockholm", "malm", "tel", "aviv", "triangeln", "durham", "the", "nya", "och", "eslov", "via", "brugge", "goteborg", "firenze", "molndal", "oslo", "ltd", "com", "www", "vaxjo", "ab"]

goodword(w) = length(w) > 1 && w ∉ badwords

function break2labels(Str)
    str = lowercase(Str)
    str = replace(str, r"ö", 'o')
    str = replace(str, r"ä", 'a')
    str = replace(str, r"å", 'a')
    str = replace(str, r"[^a-z0-9]", ' ')
    txts::Vector{String} = split(str)
    utxts = unique(strip.(txts))
    # filter!(goodword, utxts)
end

regexit(words) = isempty(words) ? r"" : Regex(string("[", join(words, "|"), "]"))

a = main()
d = a[:description]

uds = unique(break2labels(u) for u in d)

for i in d
    s = 0
    words = break2labels(i)
    for u in uds
        if
            s += 1
        end
    end
    println(s)
end

find(uds,
# a[:label] = join.(break2labels.(a[:description]), " ")

a = main()
d = a[:description]
uds = unique(break2labels(u) for u in d)
# sum(length, uds)
words = cat(1, uds...)
x = countmap(words)
w = collect(keys(x))
c = collect(values(x))
o = sortperm(c, rev=true)
y = OrderedDict(w[i] => c[i] for i in o)

x = SortedDict(x)

y = sort(x, by=last)

i = 105
d[i]
break2labels(d[i])

for j in 1:length(uds)
    u = uds[j]
    n = length(u)
    if n > 1
        for i in 1:n - 1
            if u[1:i] ∉ uds
                uds[j] = u[1:i]
            end
        end
    end
end

sum(length, uds)

a = main()

msg = """
Liability
    1) Transfer
    2) Other

Expense
    Food
        3) Cafe
        4) Groceries
        5) Restaurant

    Alcohol
        6) Store
        7) Bar

    Health
        8) Sport
        9) Other

    Transport
        10) Private
        11) Public

    FixedCost
        12) Rent
        13) Internet
        14) Electricity
        15) Phone
        16) Daycare
        17) CSN

    18) Shopping

    19) Fun

    20) Travel

    21) Other

Revenue
    22) ParentalSupport
    23) Salary
    24) ExternalTransfer
    25) Other

26) Undefined
"""

categories = [:Transfer, :Cafe, :Groceries, :Restaurant, :Systembolaget, :Bar, :Sport, :OtherHealth, :PrivateTransportation, :PublicTransportation, :Rent, :Internet, :Electricity, :Phone, :Daycare, :Electricity, :Phone, :PeriodsKort, :CSN, :Shopping, :Fun, :Travel, :OtherExpenses, :Föräldrarstöd, :Salary, :ExternalTransfer, :Undefined]
c = zeros(Int, nrow(a))
for i in 1:1#nrow(a)
    println(a[i,:])
    println(msg)
    x = readline()
    c[i] = parse(Int, x)
end

b = filter(x -> x[:category] <: Undefined, a)
open("tmp.txt", "w") do o
    println.(o, b[:description])
end

const days = Dates.Day.([0, -1, 1])
function findtransfer(a, i)::Union{Missing, Int}
    for d in days
        date = a[i, :date] + d
        i1 = findfirst(x -> x == date, a[:date])
        i1 == 0 && continue
        i2 = findlast(x -> x == date, a[:date])
        i2 == 0 && continue
        for ii in i1:i2
            if a[ii, :amount] == -a[i, :amount]
                return ii
            end
        end
    end
    return missing
end

tokill = Vector{Pair{Int, Int}}()
for i in 1:nrow(a)
    if a[i, :amount] > 0
        ii = findtransfer(a, i)
        if !ismissing(ii)
            push!(tokill, Pair(i, ii))
        end
    end
end

function findtransfer(r, a)::Union{Missing, Int}
    r[:amount] ≥ 0 && return missing
    i1 = findfirst(x -> x[:date] ≥ r[:date] - d, eachrow(a))
    i2 = findlast(x -> x[:date] ≤ r[:date] + d, eachrow(a))
    candidates = filter(i -> a[i, :amount] == -r[:amount], i1:i2)
    isempty(candidates) && return missing
    _, i = findmax(compare(Hamming(), r[:description], a[candidate, :description]) for candidate in candidates)
    candidates[i]
end

transfer(x) = r"(?:Insättning|plutt|SUGARMAMA|yakirhuvud2yakirbuffert|ninnalön2ninnahuvud|buffert|Övf|unknown|Överföring|swish|Insättning)"i(x)

i1 = findfirst(x -> x[:amount] > 0, eachrow(a))
i2 =
struct Transfer
    i1::Int
    i2::Int
    d::Dates.Day
end
ps = Dict{Float64, Vector{Transfer}}()
for i in i1:i2 - 1
    if a[i, :amount] > 0 && transfer((a[i, :description]))
        for j in i + 1:i2
            if a[i, :amount] == -a[j, :amount] && transfer(a[j, :description])
                Δ = abs(a[i, :date] - a[j, :date])
                if haskey(ps, a[j, :amount])
                    push!(ps[a[i, :amount]], Transfer(i, j, Δ))
                else
                    ps[a[i, :amount]] = [Transfer(i, j, Δ)]
                end
            end
        end
    end
end

for (k,v) in ps
    while !isempty(v)
        _, i = findmin(x.d for x in v)
        m = v[i]
        push!(tokill, Pair(m.i1, m.i2))
        filter!(x -> x.i1 ≠ m.i1 && x.i2 ≠ m.i2, v)
    end
end

tokill = Vector{Pair{Int, Int}}()
for i1 = 1:nrow(a)
    _i2 = findlast(j -> a[j, :date] ≤ a[i, :date] + d, rng)

    function distance(a, i, j)
        if a[i, :amount] == -a[j, :amount] && transfer(a[j, :description])
            abs(a[i, :date] - a[j, :date])
        else
            Dates.Day(10000)
        end
    end

    d = Dates.Day(3)
    tokill = Vector{Pair{Int, Int}}()
    for i in 1:nrow(a) - 1
        if a[i, :amount] < 0 && transfer((a[i, :description]))
            rng = i+1:nrow(a)
            _i2 = findlast(j -> a[j, :date] ≤ a[i, :date] + d, rng)
            if _i2 ≠ 0
                rng = rng[1:_i2]
                v, j = findmin(distance(a, i, j) for j in rng)
                if v ≠ Dates.Day(10000)
                    push!(tokill, i => rng[j])
                end
            end
        end
    end

    for i in 1:nrow(a)
        j
        d = [distance(a, i, j) for i = 1:nrow(a), j = 1:nrow(a)]

        for r1 in eachrow(a), r2 in eachrow(a)

            d = Dates.Day(1)
            tokill = Vector{Pair{Int, Int}}()
            for (i, r) in enumerate(eachrow(a))
                j = findtransfer(r, a)
                if !ismissing(j) && j ∉ last.(
                                              push!(tokill, i => j)
                                          end
                                      end

                                      open("tmp.txt", "w") do o
                                          for k in tokill
                                              i,j = k
                                              println(o, i, j, a[[i,j], :])
                                          end
                                      end

a = ["Comment" "Category";
     "2018-06-10 ICA SUPERMARKET SODER, MALMO" "Expense:Food:Groceries";
     "2018-06-05 SKANETRAFIKEN, MALMO C" "Expense:Transport:Public"]
=#


text:
# Abstract expression argument
const _ArgExpr{C, V} = Union{JuMP.GenericAffExpr{C, V}, JuMP.GenericQuadExpr{C, V}, V}

## Return the promoted datastructure for the indicated variable type
# GenericAffExpr
function _arg_promote(type::Type{<:GeneralVariableRef},
                      expr::JuMP.GenericAffExpr{C, V}) where {C,
                      V <: GeneralVariableRef}
    return convert(JuMP.GenericAffExpr{C, type}, expr)
end

# GenericQuadExpr
function _arg_promote(type::Type{<:GeneralVariableRef},
                      expr::JuMP.GenericQuadExpr{C, V}) where {C,
                      V <: GeneralVariableRef}
    return convert(JuMP.GenericQuadExpr{C, type}, expr)
end

# Fallback (variables, constants, etc.)
function _arg_promote(type::Type{<:GeneralVariableRef}, var::GeneralVariableRef)
    return JuMP.GenericAffExpr(0.0, Pair{type, Float64}(var, 1.0))
end

## Extend MutableArithmetics.mutable_operate to handle mixed variable types
# 2 argument adding
function _MA.mutable_operate!(::typeof(+),
                              expr::JuMP._GenericAffOrQuadExpr{C, V},
                              x::_ArgExpr{C, W}) where {C,
                              V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    expr = _arg_promote(type, expr)
    x = _arg_promote(type, x)
    return JuMP.add_to_expression!(expr, x)
end

# 2 argument adding/multiplying
function _MA.mutable_operate!(::typeof(_MA.add_mul),
                              expr::JuMP._GenericAffOrQuadExpr{C, V},
                              x::_ArgExpr{C, W}) where {C,
                              V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    expr = _arg_promote(type, expr)
    x = _arg_promote(type, x)
    return JuMP.add_to_expression!(expr, x)
end

# 2 argument substracting
function _MA.mutable_operate!(::typeof(-),
                              expr::JuMP._GenericAffOrQuadExpr{C, V},
                              x::_ArgExpr{C, W}) where {C,
                              V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    expr = _arg_promote(type, expr)
    x = _arg_promote(type, x)
    return JuMP.add_to_expression!(expr, -1, x)
end

# 2 argument substracting/multiplying
function _MA.mutable_operate!(::typeof(_MA.sub_mul),
                              expr::JuMP._GenericAffOrQuadExpr{C, V},
                              x::_ArgExpr{C, W}) where {C,
                              V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    expr = _arg_promote(type, expr)
    x = _arg_promote(type, x)
    return JuMP.add_to_expression!(expr, -1, x)
end

# 3 argument adding/multiplying (no constants)
function _MA.mutable_operate!(::typeof(_MA.add_mul),
                              expr::JuMP._GenericAffOrQuadExpr{C, V},
                              x::_ArgExpr{C, W},
                              y::_ArgExpr{C, Q}) where {C,
                              V <: GeneralVariableRef, W <: GeneralVariableRef,
                              Q <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    type = _var_type_parser(type, Q)
    expr = _arg_promote(type, expr)
    x = _arg_promote(type, x)
    y = _arg_promote(type, y)
    return JuMP.add_to_expression!(expr, x, y)
end

# 3 argument adding/multiplying (2nd constant)
function _MA.mutable_operate!(::typeof(_MA.add_mul),
                              expr::JuMP._GenericAffOrQuadExpr{C, V},
                              x::_ArgExpr{C, W},
                              y::JuMP._Constant) where {C,
                              V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    expr = _arg_promote(type, expr)
    x = _arg_promote(type, x)
    return JuMP.add_to_expression!(expr, x, y)
end

# 3 argument adding/multiplying (1st constant)
function _MA.mutable_operate!(::typeof(_MA.add_mul),
                              expr::JuMP._GenericAffOrQuadExpr{C, V},
                              x::JuMP._Constant,
                              y::_ArgExpr{C, W}) where {C,
                              V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    expr = _arg_promote(type, expr)
    y = _arg_promote(type, y)
    return JuMP.add_to_expression!(expr, x, y)
end

# 3 argument substracting/multiplying (no constants)
function _MA.mutable_operate!(::typeof(_MA.sub_mul),
                              expr::JuMP._GenericAffOrQuadExpr{C, V},
                              x::_ArgExpr{C, W},
                              y::_ArgExpr{C, Q}) where {C,
                              V <: GeneralVariableRef, W <: GeneralVariableRef,
                              Q <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    type = _var_type_parser(type, Q)
    expr = _arg_promote(type, expr)
    x = _arg_promote(type, x)
    y = _arg_promote(type, y)
    return JuMP.add_to_expression!(expr, -x, y)
end

# 3 argument substracting/multiplying (1st constant)
function _MA.mutable_operate!(::typeof(_MA.sub_mul),
                              expr::JuMP._GenericAffOrQuadExpr{C, V},
                              x::JuMP._Constant,
                              y::_ArgExpr{C, W}) where {C,
                              V <: GeneralVariableRef, W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    expr = _arg_promote(type, expr)
    y = _arg_promote(type, y)
    return JuMP.add_to_expression!(expr, -x, y)
end

# 3 argument substracting/multiplying (terminal constant)
function _MA.mutable_operate!(::typeof(_MA.sub_mul),
                              expr::JuMP._GenericAffOrQuadExpr{C, V},
                              x::_ArgExpr{C, W},
                              y::JuMP._Constant) where {C, V <: GeneralVariableRef,
                                                        W <: GeneralVariableRef}
    type = _var_type_parser(V, W)
    expr = _arg_promote(type, expr)
    x = _arg_promote(type, x)
    return JuMP.add_to_expression!(expr, x, -y)
end


text:
module MultiArmedBanditEnv

using ..Reinforce
using Distributions

export MultiArmedBandit

"""
The multi-armed bandit environment.

    MultiArmedBandit(k, n = 100; σ = 1)

- `k` is the number of available arms.
  The reward distribution for each arms is a Normal distribution centered
  at the range of [-1, 1] and the standard deviation is 1.
- `n` is the max steps for each episode.
- `σ` controls the the standard deviation for all Normal distribution.

    MultiArmedBandit(x::Vector{<:Distribution}...)

In case that you want to other distributions as the reward distribution.
"""
mutable struct MultiArmedBandit{K,D<:Vector} <: AbstractEnvironment
  arms::D
  n::Int  # max steps
  r::Float64
end

function MultiArmedBandit(k::Int, n::Int = 1000; σ::Real = 1)
  k ≤ 0 && throw(ArgumentError("k must > 0"))
  arms = map(i -> Normal(rand(Uniform(-1, 1)), σ), 1:k)
  MultiArmedBandit{k,typeof(arms)}(arms, n, 0)
end

function MultiArmedBandit(x::Vararg{Distribution,N}) where {N}
  y = collect(x)
  MultiArmedBandit{N,typeof(y)}(y, N, 0)
end

Reinforce.state(::MultiArmedBandit) = nothing
Reinforce.reward(env::MultiArmedBandit) = env.r
Reinforce.reset!(env::MultiArmedBandit) = (env.r = 0; env)
Reinforce.actions(::MultiArmedBandit{K}, s) where {K} = Base.OneTo(K)
Reinforce.step!(env::MultiArmedBandit, s, a::Int) = (env.r = rand(env.arms[a]); (env.r, nothing))
Reinforce.maxsteps(env::MultiArmedBandit) = env.n
Reinforce.ismdp(::MultiArmedBandit) = true

end  # module MultiArmedBanditEnv


text:
"""
    binary_search(array::Array{Int64,1}, target::Int64)

The binary search algorithm (BSA) finds a target value's position within a sorted array by
using a half-interval cut per each cycle. Thus, the BSA compares the target value to the
value of the array's middle element. In the case of inequality, the half array-piece in
which the target cannot be will be erased. Next, the search continues on the remaining
half array-piece and starts taking the middle element to compare it to the target value.
This procedure has to be continued until the target value is found. The search may have to
be stopped with a remaining empty half array-piece; consequently, the target is not in the
array. For more information see: [https://en.wikipedia.org/wiki/Binary_search_algorithm](https://en.wikipedia.org/wiki/Binary_search_algorithm)

# Arguments
- `array::Array{Int64,1}`: Sorted array of integers
- `target::Int64`: Target-value to find the position

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [10,11, 12, 14, 16, 18, 19, 20, 21, 22, 23, 24, 33, 35, 42, 47]
julia> target = 12
julia> ClassicAlgorithmsCollections.binary_search(arr, target)
3
```
"""
function binary_search(array::Array{Int64,1}, target::Int64)
    left = 1
    right = length(array)

    while left <= right
        # Generating the middle element position
        mid = fld((left + right), 2)
        # If element is larger than mid, then it  can only be present in right subarray
        if array[mid] < target
            left = mid + 1
            # If element is smaller than mid, then it can only be present in left subarray
        elseif array[mid] > target
            right = mid - 1
            # If element is present at the middle itself
        else
            return mid
        end
    end
    return nothing
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
# 1-byte unsigned integer -> UInt64
function size1unpack(bytes::AbstractArray{UInt8,1})
    size1unpack(view(bytes, 1:length(bytes)))
end

function size1unpack(bytes::SubArray{UInt8})
    size1unpack(bytes[1])
end

function size1unpack(byte::UInt8)
    UInt64(byte)
end

# 2-bytes unsigned integer in big endian format -> UInt64
function size2unpack(bytes::AbstractArray{UInt8,1})
    size2unpack(view(bytes, 1:length(bytes)))
end

function size2unpack(bytes::SubArray{UInt8})
    size = reinterpret(UInt16, reverse(view(bytes, max(1, length(bytes) - 2):length(bytes))))
    if length(size) > 0
        UInt64(size[1])
    else
        UInt64(0)
    end
end

# 4-bytes unsigned integer in big endian format -> UInt64
function size4unpack(bytes::AbstractArray{UInt8,1})
    size4unpack(view(bytes, 1:length(bytes)))
end

function size4unpack(bytes::SubArray{UInt8})
    size = reinterpret(UInt32, reverse(view(bytes, max(1, length(bytes) - 4):length(bytes))))
    if length(size) > 0
        UInt64(size[1])
    else
        UInt64(0)
    end
end

function int1unpack(bytes::AbstractArray{UInt8,1})
    int1unpack(view(bytes, 1:length(bytes)))
end

function int1unpack(bytes::SubArray{UInt8})
    int1unpack(bytes[1])
end

function int1unpack(byte::UInt8)
    Int(byte)
end

function int4unpack(bytes::AbstractArray{UInt8,1})
    int4unpack(view(bytes, 1:length(bytes)))
end

function int4unpack(bytes::SubArray{UInt8})
    Int(reinterpret(Int32, reverse(view(bytes, max(1, length(bytes) - 4):length(bytes))))[1])
end

function floatunpack(bytes::AbstractArray{UInt8,1})
    floatunpack(view(bytes, 1:length(bytes)))
end

function floatunpack(bytes::SubArray{UInt8})
    reinterpret(Float64, reverse(view(bytes, max(1, length(bytes) - 8):length(bytes))))[1]
end

function charintpack(value::Integer, size::Int)
    bytes = zeros(UInt8, size)
    for i in 1:size
        bytes[i] = UInt8(value & 0xff)
        value = value >>> 8
    end
    reverse(bytes)
end

function charintpack(value::Integer)
    bytes = []
    while value != 0
        bytes = vcat(bytes, UInt8(value & 0xff))
        value = value >>> 8
    end
    bytes
end

function charint4pack(integer::Integer)
    charintpack(integer, 4)
end

function charint2pack(integer::Integer)
    charintpack(integer, 2)
end

function charsignedint4pack(integer::Integer)
    charintpack(integer, 4)
end

function lencheck(bytes::AbstractArray{UInt8,1}, limit::Int)
    lencheck(view(bytes, 1:length(bytes)), limit)
end

function lencheck(bytes::SubArray{UInt8}, limit::Int)
    len = UInt64(length(bytes))
    lencheck(len, len < limit, bytes)
end

function lencheck(len::UInt64, limit::UInt64, bytes::AbstractArray{UInt8,1})
    lencheck(len, limit, view(bytes, 1:length(bytes)))
end

function lencheck(len::UInt64, limit::UInt64, bytes::SubArray{UInt8})
    lencheck(limit, len < limit, bytes)
end

function lencheck(len::UInt64, pred::Bool, bytes::AbstractArray{UInt8,1})
    lencheck(len, pred, view(bytes, 1:length(bytes)))
end

function lencheck(len::UInt64, pred::Bool, bytes::SubArray{UInt8})
    if pred
        throw(IncompleteData(bytes))
    end
    len
end

function compressterm(encodedterm, compression::Bool)
    compressterm(encodedterm, 6)
end

function compressterm(encodedterm, compression::Int)
    throw(NotImplemented())
end

function decompressterm(bytes::AbstractArray{UInt8,1})
    throw(NotImplemented())
end