FEM_frame.jl

# Author:  JV Carstensen, CEE, MIT (JK Guest, Civil Eng, JHU)
# Revised: Aug 22 2017, JVC
# Revised: Amira Abdel-Rahman
# =======================================================================

function FEM_frame(data,A)

    # FEM solver for frame elements
    # Author:  JV Carstensen, CEE, MIT (JK Guest, Civil Eng, JHU)
    # Revised: Aug 22 2017, JVC
    # Revised: Amira Abdel-Rahman

    # =======================================================================



    # ---- READ IN DATA -----------------------------------------------------
    # read in external file
    E,f,dis,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te= data;
    g=zeros(size(dis))
    g.=dis


    # ---- NUMBER THE EQUATIONS ---------------------------------------------
    # line 380
    id,neq = number_eq(idb,ndf,nnp);


    # ---- FORM THE ELEMENT STIFFNESS MATRICES ------------------------------
    # line 324
    nee = ndf*nen;                   # number of element equations
    Ke  = zeros(nee,nee,nel);
    # Te  = zeros(nen*1,nen*nsd,nel);  # *1 is specific to frame
    Te  = zeros(nee,nee,nel);
    for i = 1:nel
        Ke[:,:,i],Te[:,:,i] = Ke_frame(A[i],E[i],ien[:,i],nee,nsd,xn);
    end

    # display("Element Stiffness Matrix")
    # display(Ke)


    # ---- PERFORM GLOBAL TO LOCAL MAPPING ----------------------------------
    # line 237
    LM  = zeros(nee,nel);
    for i = 1:nel
        LM[:,i] = get_local_id(id,ien[:,i],ndf,nee,nen);
    end


    # ---- IF THERE IS FREE DEGREES OF FREEDOM, THEN SOLVE THE EQUILIBRIUM --
    if (neq > 0)

        # get global force vector (line 275)
        F = globalF(f,g,id,ien,Ke,LM,ndf,nee,nel,nen,neq,nnp);

        # solve Kd - F = 0 (line 521)
        d,K = solveEQ(F,LM,Ke,nee,nel,neq);
    end

    # display("F")
    # display(F)


    # ---- POST-PROCESS THE RESULTS -----------------------------------------
    # line 419
    dcomp,axial,stress,strain,Fe = post_processing(A,d,E,g,id,ien,Ke,ndf,nee,nel,nen,nnp,Te);
    # display("dcomp")
    # display(dcomp)

    # display("axial")
    # display(axial)

    # display("stress")
    # display(stress)
    # ---- COMPUTE REACTION FORCES ------------------------------------------
    # line 482
    Rcomp,idr = reactions(idb,ien,Fe,ndf,nee,nel,nen,nnp);
    # display("Rcomp")
    # display(Rcomp)

    # ---- PLOT THE STRUCTURE -----------------------------------------------
    # read in external file

    # set the plot factor for the thickness of the frame elements
    #  plot_fac_bar = 1/A[i];
    #  A_min        = 0;
    #  plot_frame(A,A_min,f,idr,ien,nel,nnp,nsd,plot_fac_bar,xn);
    return K,F,d,stress,dcomp,g
    # =======================================================================
        
end
# =======================================================================
function add_d2dcomp(dcomp,d,id,ndf,nnp)

    # function that adds the displacements of the free degrees of freedom to
    # the nodal displacements
    # -----------------------------------------------------------------------
    # dcomp(ndf,nnp)  = nodal displacements
    # d(neq,1)        = displacement at free degrees of freedom
    # id(ndf,nnp)     = equation numbers of degrees of freedom
    # ndf             = number of degrees of freedom per node
    # nnp             = number of nodal points

    #------------------------------------------------------------------------

    # loop over nodes and degrees of freedom
    for n=1:nnp
        for i=1:ndf

            # if it is a free dof then add the global displacement
            if (id[i,n]>0)       
                dcomp[Int(i),Int(n)] = dcomp[Int(i),Int(n)]+d[Int(id[Int(i),Int(n)])];
            end
        end
    end
    return dcomp
end
# =======================================================================

function add_loads_to_force(F,f,id,ndf,nnp)

    # function that adds nodal forces to the global force vector
    # -----------------------------------------------------------------------
    # F(neq,1)        = global force vector
    # f(ndf,nnp)      = prescribed nodal forces
    # id(ndf,nnp)     = equation numbers of degrees of freedom
    # ndf             = number of degrees of freedom per node
    # nnp             = number of nodal points

    #------------------------------------------------------------------------

    # loop over nodes and degrees of freedom
    for n = 1:nnp
        for i = 1:ndf

            # get the global equation number 
            M = id[i,n];  

            # if free degree of freedom, then add nodal load to global force 
            # vector
            if (M > 0)        
                F[Int(floor(M))] = F[Int(floor(M))] + f[Int(floor(i)),Int(floor(n))];
            end
        end
    end
    return F

end
# =======================================================================

function addforce(F,Fe,LM,nee)

    # function that adds element forces to the global force vector
    # -----------------------------------------------------------------------
    # F(neq,1)        = global force vector
    # Fe(nee,1)       = element force vector
    # LM(nee,nel)     = global to local map for the element
    # nee             = number of element equations

    # =======================================================================

    # loop over rows of Fe
    for i = 1:nee

        # get the global equation number for local equation i
        M = LM[i];     

        # if free dof (eqn number > 0) add to F vector
        if (M > 0)         
            F[Int(M)]=F[Int(M)]+Fe[Int(i)];
        end
    end
    return F


end
# =======================================================================

function get_de_from_dcomp(dcomp,ien,ndf,nen)

    #  extracts element displacement vector from complete displacement vector
    #------------------------------------------------------------------------
    # dcomp(ndf,nnp)  = nodal displacements
    # ien(nen,1)      = element connectivity
    # ndf             = number of degrees of freedom per node
    # nen             = number of element equations
    #
    # de(nen,1)       = element displacements
    # 
    #------------------------------------------------------------------------
    
    de = zeros((nen-1)*ndf+ndf,1);

    # loop over number of element nodes
    for i = 1:nen

        # loop over number of degrees of freedom per node
        for j = 1:ndf

            # get the local element number and place displacement in de
            leq     = (i-1)*ndf+j;  
            
            de[leq] = dcomp[j,Int(floor(ien[i]))];
        end
    end
    return de
end
# =======================================================================

function get_local_id(id,ien,ndf,nee,nen)

    # functions that performs the global to local mapping of equation numbers
    # -----------------------------------------------------------------------
    # id(ndf,nnp)     = equation numbers of degrees of freedom
    # ien(nen,1)      = element connectivity
    # ndf             = number of degrees of freedom per node
    # nee             = number of element equations
    # nen             = number of element equations
    #
    # LM(nee,1)       = global to local map for element

    # =======================================================================

    # initialize global-local mapping matrix
    LM  = zeros(nee,1);        

    # initialize local equation number counter
    k = 0;                      

    # loop over element nodes
    for i = 1:nen         

        # loop over degrees of freedom at each node
        for j = 1:ndf

            # update counter and prescribe global equation number
            k     = k+1;
            LM[k] = id[j,Int(floor(ien[i]))];
        end
    end
    return LM
end
# =======================================================================

function globalF(f,g,id,ien,Ke,LM,ndf,nee,nel,nen,neq,nnp) 

    # function that assembles the global load vector
    # -----------------------------------------------------------------------
    # id(ndf,nnp)     = equation numbers of degrees of freedom
    # f(ndf,nnp)      = prescribed nodal forces
    # g(ndf,nnp)      = prescribed nodal displacements
    # ien(nen,nel)    = element connectivities
    # Ke(nee,nee,nel) = element stiffness matrices
    # ndf             = number of degrees of freedom per node
    # nee             = number of element equations
    # nel             = number of elements
    # nen             = number of element equations
    # neq             = number of equations
    # nnp             = number of nodal points
    #
    # F(neq,1)        = global force vector

    # =======================================================================

    # initialize
    F = zeros(neq,1);

    # Insert applied loads into F
    F = add_loads_to_force(F,f,id,ndf,nnp);

    # Compute forces from applied displacements (ds~=0) and insert into F
    Fse = zeros(nee,nel);  

    # loop over elements
    for i = 1:nel

        # get dse for current element
        dse      = get_de_from_dcomp(g,ien[:,i],ndf,nen);  

        # compute element force
        Fse[:,i] = -Ke[:,:,i]*dse;

        # assemble elem force into global force vector
        F        = addforce(F,Fse[:,i],LM[:,i],nee);     
    end

    return F
end
# =======================================================================

function  Ke_frame(A,E,ien,nee,nsd,xn)

    # function that computes the global element stiffness matrix for a frame
    # element
    # -----------------------------------------------------------------------
    # A(1,1)          = cross-sectional area of elements
    # E(1,1)          = Young's modulus of elements
    # ien(nen,1)      = element connectivity
    # nee             = number of element equations
    # nsd             = number of spacial dimensions
    # xn(nsd,nnp)     = nodal coordinates
    #
    # Ke(nee,nee,1)   = global element stiffness matrix
    # Te(nee,nee,1)   = global to local transformation matrix for element

    # =======================================================================

    # form vector along axis of element using nodal coordinates
    v = xn[:,Int(floor(ien[2]))]-xn[:,Int(floor(ien[1]))];
    

    # compute the length of the element
    Le = norm(v,2);

    # rotation of parent domain
    #   rot=[ cos(theta_x)  cos(theta_y)  cos(theta_z) ]'
    rot = v/Le;

    l=rot[1]
    m=rot[2]
    n=rot[3]

    D=sqrt(l^2+m^2+n^2)
    # D=1.0

    b=sqrt(A)
    h=sqrt(A)

    G=E * 1 / 2 #todo shear_modulus
    ixx = b*h^3/12
    iyy = b*h^3/12
    j=b*h*(b*b+h*h)/12;#todo check
    l0=Le
    l02 = l0 * l0
    l03 = l0 * l0 * l0
    

    # local element stiffness matrix
    # ke = E*A/Le*[  1  -1
    #     -1   1 ];

    ke = [  [E*A/l0  0  0  0  0  0  -E*A/l0  0  0  0  0  0];
            [0  12*E*ixx/l03  0  0  0  6*E*ixx/l02  0  -12*E*ixx/l03  0  0  0  6*E*ixx/l02];
            [0  0  12*E*iyy/l03  0  -6*E*iyy/l02  0  0  0  -12*E*iyy/l03  0  -6*E*iyy/l02  0];
            [0  0  0  G*j/l0  0  0  0  0  0  -G*j/l0  0  0];
            [0  0  -6*E*iyy/l02  0  4*E*iyy/l0  0  0  0  6*E*iyy/l02  0  2*E*iyy/l0  0];
            [0  6*E*ixx/l02  0  0  0  4*E*ixx/l0  0  -6*E*ixx/l02  0  0  0  2*E*ixx/l0];
            [-E*A/l0  0  0  0  0  0  E*A/l0  0  0  0  0  0];
            [0  -12*E*ixx/l03  0  0  0  -6*E*ixx/l02  0  12*E*ixx/l03  0  0  0  -6*E*ixx/l02];
            [0  0  -12*E*iyy/l03  0  6*E*iyy/l02  0  0  0  12*E*iyy/l03  0  6*E*iyy/l02  0];
            [0  0  0  -G*j/l0  0  0  0  0  0  G*j/l0  0  0];
            [0  0  -6*E*iyy/l02  0  2*E*iyy/l0  0  0  0  6*E*iyy/l02  0  4*E*iyy/l0  0];
            [0  6*E*ixx/l02  0  0  0  2*E*ixx/l0  0  -6*E*ixx/l02  0  0  0  4*E*ixx/l0]];


    # Transformation matrix: global to local coordinate system
    if (nsd == 2)   # 2D case

        # # frame Te is nen x ndf*nen array
        # Te = [ rot[1]  rot[2]       0       0
        #             0       0  rot[1]  rot[2] ];
        # Frame Te is nen x ndf*nen array
        Te = [ rot[1]  rot[2]   0          0       0  0
              -rot[2]  rot[1]   0          0       0  0
                   0       0    1          0       0  0
                   0       0    0     rot[1]  rot[2]  0
                   0       0    0    -rot[2]  rot[1]  0
                   0       0    0          0       0  1];

    elseif (nsd == 3)   # 3D case

        # # frame Te is nen x ndf*nen array
        # Te = [ rot[1]  rot[2]  rot[3]       0       0       0  
        #             0       0       0  rot[1]  rot[2]  rot[3] ];
        # Frame Te is nen x ndf*nen array
        # Te = [ rot[1]  rot[2]  rot[3]       0       0       0  
        #             0       0       0  rot[1]  rot[2]  rot[3] ];


        Te = [       l       m  n      0       0 0       0       0 0       0      0 0
                  -m/D     l/D  0      0       0 0       0       0 0       0      0 0
                -l*n/D  -m*n/D  D      0       0 0       0       0 0       0      0 0
                     0       0  0      l       m n       0       0 0       0      0 0
                     0       0  0   -m/D     l/D 0       0       0 0       0      0 0
                     0       0  0 -l*n/D  -m*n/D D       0       0 0       0      0 0
                     0       0  0      0       0 0       l       m n       0      0 0
                     0       0  0      0       0 0    -m/D     l/D 0       0      0 0
                     0       0  0      0       0 0  -l*n/D  -m*n/D D       0      0 0
                     0       0  0      0       0 0       0       0 0       l      m n
                     0       0  0      0       0 0       0       0 0    -m/D    l/D 0
                     0       0  0      0       0 0       0       0 0  -l*n/D -m*n/D D];
    end
    # compute the global element stiffness matrix
    Ke = zeros(nee,nee);
    Ke = Te'*ke*Te;
    # println(size(Ke))
    # println(size(Te))

    return Ke,Te
end
# =======================================================================

function number_eq(idb,ndf,nnp)

    # function that numbers the unknown degrees of freedom (equations)
    # -----------------------------------------------------------------------
    # idb(ndf,nnp)    = 1 if the degree of freedom is prescribed, 0 otherwise
    # ndf             = number of degrees of freedom per node
    # nnp             = number of nodal points
    #
    # id(ndf,nnp)     = equation numbers of degrees of freedom
    # neq             = number of equations (tot number of degrees of freedom)

    # =======================================================================

    # initialize id and neq
    id  = zeros(ndf,nnp);  
    neq = 0;              

    # loop over nodes
    for n = 1:nnp

        # loop over degrees of freedom
        for i = 1:ndf
            if idb[i,n] == 0 

                # udate # of equations
                neq = neq + 1;      

                # if no prescribed displacement at dof i of node n
                #   => give an equation # different from 0
                id[i,n] = neq;      

            end
        end
    end
    return id,neq
end

# =======================================================================
function post_processing(A,d,E,g,id,ien,Ke,ndf,nee,nel,nen,nnp,Te)

    # function that performs post processing for frame elements
    # -----------------------------------------------------------------------
    # A(nel,1)        = cross-sectional area of elements
    # d(neq,1)        = displacements at free degrees of freedom
    # E(nel,1)        = Young's modulus of elements
    # g(ndf,nnp)      = prescribed nodal displacements
    # id(ndf,nnp)     = equation numbers of degrees of freedom
    # ien(nen,nel)    = element connectivities
    # Ke(nee,nee,nel) = element stiffness matrices
    # ndf             = number of degrees of freedom per node
    # nee             = number of element equations
    # nel             = number of elements
    # nen             = number of element equations
    # nnp             = number of nodes
    # Te(nee,nee,nel) = element transformation matrices
    #
    # dcomp(ndf,nnp)  = nodal displacements
    # axial(nel,1)    = axial element forces
    # stress(nel,1)   = element stresses
    # strain(nel,1)   = element strains
    # Fe(nee,nel)     = element forces

    # =======================================================================

    # get the total displacement of the structure in matrix form dcomp(nsd,nnp)
    dcomp = add_d2dcomp(g,d,id,ndf,nnp);

    # initalize evaluation of global element forces Fe, local element forces 
    # fe, axial forces, element stresses and strains
    Fe     = zeros(nee,nel);      
    de     = zeros(nee,nel);  
    fe     = zeros(nee,nel);  # element local force vector 
    axial  = zeros(nel,1); 
    stress = zeros(nel,1); 
    strain = zeros(nel,1); # element axial, stress, strain

    

    # loop over elements
    for i=1:nel

        # get the element displacaments
        de[:,i] = get_de_from_dcomp(dcomp,ien[:,i],ndf,nen);

        # compute the element forces
        Fe[:,i] = Ke[:,:,i]*de[:,i];

        # transform Fe to the local coordinate system
        fe[:,i] = Te[:,:,i]*Fe[:,i];

        # Compute the axial force, stress, strain
        axial[i] = fe[7,i] ;         # Use second entry for frame element
        stress[i] = axial[i]/A[i];   
        strain[i] = stress[i]/E[i];
    end

    # display(fe[:,1])

    # display("Fe")
    # display(Fe)
    return dcomp,axial,stress,strain,Fe
end
# =======================================================================

function reactions(idb,ien,Fe,ndf,nee,nel,nen,nnp)

    # function that computes the reaction forces on the structure
    # -----------------------------------------------------------------------
    # idb(ndf,nnp)    = 1 if the degree of freedom is prescribed, 0 otherwise
    # ien(nen,nel)    = element connectivities
    # Fe(nee,nel)     = element forces
    # ndf             = number of degrees of freedom per node
    # nee             = number of element equations
    # nel             = number of elements
    # nen             = number of element equations
    # nnp             = number of nodes
    #
    # Rcomp(ndf,nnp)  = nodal reactions
    # idbr(ndf,nnp)   = 0 if the degree of freedom is prescribed,  otherwise

    # =======================================================================

    # switch BC marker and number the equations for the reaction forces
    idbr = 1 .- idb;  
    idr,neqr = number_eq(idbr,ndf,nnp);

    # assemble reactions R from element force vectors Fe 
    R   = zeros(neqr,1);  
    for i = 1:nel
        LMR = get_local_id(idr,ien[:,i],ndf,nee,nen);
        R   = addforce(R,Fe[:,i],LMR,nee);
    end

    # organize the reactions in matrix array Rcomp(ndf,nnp)
    Rcomp = zeros(ndf,nnp);
    Rcomp = add_d2dcomp(Rcomp,R,idr,ndf,nnp);

    


    return Rcomp,idr
end
# =======================================================================

function addstiff(K,Ke,LM,nee)
    
    # function that solves the equilibrium condition
    # -----------------------------------------------------------------------
    # K(neq,neq)      = global stiffness matrix
    # Ke(nee,nee,1)   = element stiffness matrix
    # LM(nee,nel)     = global to local map for the element
    # nee             = number of element equations

    # =======================================================================

    # loop over rows of Ke
    for i =1:nee

        # loop over columns of Ke

        for j = 1:nee                        

            Mr = LM[i];
            Mc = LM[j];

            if(Mr > 0 && Mc > 0)             
                # if equation #'s are non-zero add element contribution to the 
                # stiffness matrix
                K[Int(Mr),Int(Mc)] = K[Int(Mr),Int(Mc)] +Ke[Int(i),Int(j)];
            end

        end
    
    end
    
    return K
    

end 

# =======================================================================

function solveEQ(F,LM,Ke,nee,nel,neq) 

    # function that solves the equilibrium condition
    # -----------------------------------------------------------------------
    # F(neq,1)        = global force vector
    # LM(nee,nel)     = global to local maps
    # Ke(nee,nee,nel) = element stiffness matrices
    # nee             = number of element equations
    # nel             = number of elements
    # neq             = number of equations
    #
    # d(neq,1)        = displacements at free degrees of freedom

    # =======================================================================

    # assemble global stiffness matrix
    # K = zeros(neq,neq);   # Use 'sparse' for more efficient memory usage
    K=spzeros(neq,neq)
    for i = 1:nel
        K = addstiff(K,Ke[:,:,i],LM[:,i],nee);
    end

    # display("K")
    # display(K)

    # solve the equlibrium
    d = K\F;

    # display("d")
    # display(d)

    return d,K
    
end
# =======================================================================

function mapp(value, x1, y1, x2, y2)
	return (value - x1) * (y2 - x2) / (y1 - x1) + x2;
end
# =======================================================================

function plotFrame(problem,X,scale,threshold=0)
    nel=length(X)
    E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
    K,F,d,stress,dcomp,g=FEM_frame(problem,X);
    p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
    for i in 1:nel
        if X[i]>threshold
            if threshold>0
                p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                    [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                    color=RGB(0.0,0.0,0.0),
                    linewidth = 3.0)
                    
            else
                p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                    [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                    color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
                    linewidth = X[i]*scale)
                    # linewidth = mapp(X[i], minimum(X[:]), maximum(X[:])+0.0001, scale/10.0, scale))
                    # p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
            end
        else
            p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                    [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                    color=RGB(1.0,1.0,1.0),
                    linewidth = 0.0)
        end

    end
    # plot!(axis=nothing,ticks=nothing, border=nothing)
    p
end
# =======================================================================

function getSetup(fileName)
    setup = Dict()
    open(fileName, "r") do f
        dicttxt = read(f,String)  # file information to string
        setup=JSON.parse(dicttxt)  # parse and transform data
    end
    return setup
    
end

# =======================================================================

function plotFrameDeformed(problem,X,scale,threshold=0,exageration=100.0)
    nel=length(X)
    E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
    K,F,d,stress,dcomp,g=FEM_frame(problem,X);
    p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
    xnn=zeros(size(xn))
    xnn.=xn .+ exageration.*dcomp
    for i in 1:nel
        if X[i]>threshold
            if threshold>0
                p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                    [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                    color=RGB(0.0,0.0,0.0),
                    linewidth = 3.0)
                p=plot!([xnn[1,Int(ien[1,i])],xnn[1,Int(ien[2,i])]],
                    [xnn[2,Int(ien[1,i])],xnn[2,Int(ien[2,i])]],label="",
                    color=RGB(0.0,1.0,0.0),
                    linewidth = 3.0)
                    
            else
                p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                    [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                    color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
                    linewidth = X[i]*scale)
                    # linewidth = mapp(X[i], minimum(X[:]), maximum(X[:])+0.0001, scale/10.0, scale))
                    # p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
            end
        else
            p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                    [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                    color=RGB(1.0,1.0,1.0),
                    linewidth = 0.0)
        end

    end
    # plot!(axis=nothing,ticks=nothing, border=nothing)
    p
end

# =======================================================================

function plotFrameDeformed3D(problem,X,scale,threshold=0,exageration=100.0)
    nel=length(X)
    E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
    K,F,d,stress,dcomp,g=FEM_frame(problem,X);
    p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
    xnn=zeros(size(xn))
    xnn.=xn .+ exageration.*dcomp[1:3, :]
    for i in 1:nel
        if X[i]>threshold
            if threshold>0
                p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                        [xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
                        [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                        color=RGB(0.0,0.0,0.0),
                        linewidth = 3.0)
                p=plot!([xnn[1,Int(ien[1,i])],xnn[1,Int(ien[2,i])]],
                        [xnn[3,Int(ien[1,i])],xnn[3,Int(ien[2,i])]],
                        [xnn[2,Int(ien[1,i])],xnn[2,Int(ien[2,i])]],label="",
                        color=RGB(0.0,1.0,0.0),
                        linewidth = 3.0)
                    
            else
                p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                        [xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
                        [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                        color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
                        linewidth = X[i]*scale)
                        # linewidth = mapp(X[i], minimum(X[:]), maximum(X[:])+0.0001, scale/10.0, scale))
                        # p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
            end
        else
            p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                    [xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
                    [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                    color=RGB(1.0,1.0,1.0),
                    linewidth = 0.0)
        end

    end
    # plot!(axis=nothing,ticks=nothing, border=nothing)
    p
end

# =======================================================================

function plotFrame3D(problem,X,scale,threshold=0)
    nel=length(X)
    E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
    K,F,d,stress,dcomp,g=FEM_frame(problem,X);
    p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
    for i in 1:nel
        if X[i]>threshold
            if threshold>0
                p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                        [xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
                        [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                        color=RGB(0.0,0.0,0.0),
                        linewidth = 3.0)
                    
            else
                p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                        [xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
                        [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                        color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
                        linewidth = X[i]*scale)
                        # linewidth = mapp(X[i], minimum(X[:]), maximum(X[:])+0.0001, scale/10.0, scale))
                        # p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
            end
        else
            p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
                    [xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
                    [xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
                    color=RGB(1.0,1.0,1.0),
                    linewidth = 0.0)
        end

    end
    # plot!(axis=nothing,ticks=nothing, border=nothing)
    p
end

# =======================================================================