FEM.jl

function FEM(data_FEM,ρ)

    # FEM solver for 2D elements
    # Author:  JV Carstensen, CEE, MIT (JK Guest, Civil Eng, JHU)
    # Revised: Oct 26 2017, JVC
    # Revised: Amira Abdel-Rahman
    
    # =======================================================================

    
    
    # ---- READ IN DATA -----------------------------------------------------
    # read in external file
    E,f,dis,idb,ien,iplane,ndf,nel,nen,nnp,nsd,snu,t,xn,Ke0 = data_FEM;
    g=zeros(size(dis))
    g.=dis
    
    
    # ---- NUMBER THE EQUATIONS ---------------------------------------------
    id,neq = number_eq(idb,ndf,nnp);
    
    
    # ---- FORM THE ELEMENT STIFFNESS MATRICES ------------------------------
    nee = ndf*nen;                   # number of element equations
    Ke  = zeros(nee,nee,nel);
    Imx_nsd  = Matrix(1I, nsd, nsd);
    zero_nsd = zeros(nsd,nsd); 
    for i = 1:nel
        Ke[:,:,i] .= (ρ[i]).*Ke_quad4(iplane,E,snu,t,nee,nen,nsd,ndf,ien[:,i],xn,Imx_nsd,zero_nsd);
    end
    
    
    # ---- PERFORM GLOBAL TO LOCAL MAPPING ----------------------------------
    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(Ke)

    de,Fe=post_processing(d,E,g,id,ien,Ke,ndf,nee,nel,nen,nnp)
    
    return K,F,d,de,Fe
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[i,n] = dcomp[i,n]+d[Int(id[i,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 = Int(id[i,n]);  
            
            # if free degree of freedom, then add nodal load to global force 
            # vector
            if (M > 0)        
                F[M] = F[M]+f[i,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 = Int(LM[i]);     
        
        # if free dof (eqn number > 0) add to F vector
        if (M > 0)         
            F[M]=F[M]+Fe[i];
        end
    end
    return F
end
    # =======================================================================
    
    
    # =======================================================================
    
function addstiff(K,Ke,LM,nee)
    
    # function that adds the element stiffness matrices to the global 
    # stiffness matrix
    # -----------------------------------------------------------------------
    # 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 = Int(LM[i]);
            Mc = Int(LM[j]);
            
            if (Mr > 0 && Mc > 0)             
                # if equation #'s are non-zero add element contribution to the 
                # stiffness matrix
                K[Mr,Mc] = K[Mr,Mc]+Ke[i,j];
            end
        end
    end
    return K
end
    # =======================================================================
    
    
    # =======================================================================
    
function B_2d_elastic(dNdx,dNdy,nen,ndf)
    
    #  Computes the strain-displacement matrix for an elastic 2d problem
    # -----------------------------------------------------------------------
    # B(3,nen*ndf)    = strain displacement matrix at current gauss point
    # dNdx            = derivative of shape function at current gauss point
    # dNdy            = derivative of shape function at current gauss point
    # nen             = number of nodes per element
    # ndf             = number of degrees of freedom per node
    #
    # =======================================================================
    
    # Initialize
    B = zeros(3,nen*ndf);
    
    # Compute B
    for i = 1:nen
    B[:,(i-1)*ndf+1:ndf*i] = [ dNdx[i]       0
                                    0   dNdy[i]
                               dNdy[i]  dNdx[i] ];
    end 
    
    # =======================================================================
    return B
end
    # =======================================================================
    
    
    # =======================================================================
    
function D_2d(E,snu,iplane,nstr)
    
    #  Computes the elastic constitutive matrix for a 2d problem
    # -----------------------------------------------------------------------
    # E               = Young's modulus
    # snu             = Poisson's ratio
    # iplane          = 1 - plane strain, 2 - plane stress
    # nstr            = number of independent stress components
    #
    # D(nstr,nstr)    = elastic constitutive matrix for element e
    #
    # =======================================================================
    
    # initialize and define D
    D = zeros(Int(nstr),Int(nstr));
    
    if (iplane == 2)         # plane stress
        coeff  = E/(1-(snu*snu));
        D[1,1] = coeff;
        D[2,2] = D[1,1];
        D[1,2] = coeff*snu;
        D[2,1] = D[1,2];
        D[3,3] = coeff*(1-snu)/2;
        
    elseif (iplane == 1)     # plane strain
        coeff  = E/((1+snu)*(1-2*snu));
        D[1,1] = coeff*(1-snu);
        D[2,2] = D[1,1];
        D[1,2] = coeff*snu;
        D[2,1] = D[1,2];
        D[3,3] = coeff*(1-2*snu)/2;
    end
    return D
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,ien[i]];
        end
    end
    return de
end
    # =======================================================================
    
    
    # =======================================================================
    
function  gausspoints_quad4(ndf,nsd)
    
    # function that defines integration points and weigths for quad4 elements
    #------------------------------------------------------------------------
    # ndf             = number of degrees of freedom per node
    # nsd             = number of spacial dimensions
    #
    # Nint            = number of gauss points in x and y
    # point           = location of gauss points in x and y
    # weight          = weigthing of gauss points for numerical integration
    # nstre           = number of independent stress components
    # nstr1           = total number of stress components
    # ngp             = number of gauss points per element
    # 
    #------------------------------------------------------------------------
    Nint=zeros(Int,2)
    Nint[1] = 2; # number of gauss points in x
    Nint[2] = 2; # number of gauss points in y
    
    # total number of gauss points
    ngp = Nint[1]*Nint[2];
    
    # Initialize
    point   = zeros(Nint[2],Nint[1]);
    weight =  zeros(Nint[2],Nint[1]);
    weight .= point;
    
    # get gauss point locations and weights for integration
    point[1,1] = (-0.577350269189626);
    point[1,2] = (-0.577350269189626);
    point[2,1] = ( 0.577350269189626);
    point[2,2] = ( 0.577350269189626);
    weight       .= weight .+1;
    
    # define number of individual stress/strain components
    nstre = nsd*(ndf+1)/2;
    
    # define number of stress components in each element (4 in 2d)
    nstr1 = 4;
    return Nint,point,weight,nstre,nstr1,ngp
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,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  jacobian_2d(dNdr,dNds,xn,nen,Imx,Jaco)
    
    #  computes the Jacobian, its determinate and inverse at current gauss
    #  point in element e
    # -----------------------------------------------------------------------
    # dNdr            = derivative of shape function at current gauss point
    # dNds            = derivative of shape function at current gauss point
    # xn(nsd,nen)     = nodal coordinates for element e
    # nen             = number of nodes per element
    # Imx             = eye(nsd) 
    # Jaco            = zeros(nsd)
    #
    # Jaco            = Jacobian matrix for current gauss point
    # detJ            = determinant of the Jacobian for current gauss point
    # InvJ            = inverse of the Jacobian matrix for current gauss point
    # =======================================================================
    
    # Add to Jacobian
    for i = 1:nen
        #for j = 1:nsd
            Jaco[1,1] = Jaco[1,1]+dNdr[i]*xn[1,i];
            Jaco[1,2] = Jaco[1,2]+dNdr[i]*xn[2,i];
            Jaco[2,1] = Jaco[2,1]+dNds[i]*xn[1,i];
            Jaco[2,2] = Jaco[2,2]+dNds[i]*xn[2,i];
        #end
    end
    
    # Find its determinant
    detJ = det(Jaco);
    
    # Find its inverse
    invJ = Imx/Jaco;
    
    return Jaco,detJ,invJ
end
    # =======================================================================
    
    
    # =======================================================================
    
function  ke_elem_quad4(Nint,point,weight,ien, xn,ndf,nen,nee,t,D,Imx_nsd,zero_nsd)
    
    #  loops over element gauss points and computes element stiffness matrix
    # -----------------------------------------------------------------------
    # Nint            = number of gauss points in x and y
    # point           = location of gauss points
    # weight          = weigthing of gauss points in numerical integration
    # ien(nen,nel)    = element connectivity
    # xn(nsd,nnp)     = nodal coordinates 
    # nen             = number of nodes per element
    # nee             = number of element equations
    # t               = element thickness
    # D               = elastic constititive matrix for the element
    # Imx             = eye(nsd) 
    # Jaco            = zeros(nsd)
    #
    # kee(nee,nee,1)  = element stiffness matrix for quad 4 element
    # =======================================================================
    
    # =======================================================================
    
    
    # Initialize 
    kee   = zeros(nee,nee);
    coeff = 0;
    
    
    for i = 1:Nint[1]
        ri = point[i,1];
        wi = weight[i,1];
        
        for j = 1:Nint[2]
            sj = point[j,2];
            wj = weight[j,2];
            sN,dNdx,dNdy,detJ,invJ = shape_quad4(ri,sj,xn,ien,nen,Imx_nsd,zero_nsd);
            
            B = B_2d_elastic(dNdx,dNdy,nen,ndf);
            
            coeff = t*wi*wj*detJ;
            
            # Compute kee = kee + transpose(B)*D*B*coeff
            kee .= kee .+(B')*D*B*coeff;
        end
        
    end
    return kee
end
    # ======================================================================= 
    
    
    # =======================================================================
    
function Ke_quad4(iplane,E,snu,t,nee,nen,nsd,ndf,ien,xn, Imx_nsd,zero_nsd)
    
    #  computes element stiffness matrix
    # -----------------------------------------------------------------------
    # iplane          = 1 - plane strain, 2 - plane stress 
    # E               = Young's modulus
    # snu             = Poisson's ratio
    # t               = element thickness
    # nee             = number of element equations
    # nen             = number of nodes per element
    # nsd             = number of spacial dimmensions
    # ndf             = number of degrees of freedom
    # weight          = weigthing of gauss points in numerical integration
    # ien(nen,nel)    = element connectivity
    # xn(nsd,nnp)     = nodal coordinates 
    # Imx             = eye(nsd) 
    # Jaco            = zeros(nsd)
    #
    # kee(nee,nee,1)  = element stiffness matrix for quad 4 element
    # =======================================================================
    
    # Define location and weights of integration points
    Nint,point,weight,nstre,nstr1,ngp = gausspoints_quad4(ndf,nsd);
    
    # get constitutive matrix
    nstr = nsd*(nsd+1)/2;
    D = D_2d(E,snu,iplane,nstr);
    
    
    # Compute element matrix
    kee = ke_elem_quad4(Nint,point,weight,ien,xn,ndf,nen,nee,t,D,Imx_nsd,zero_nsd);
    
    return kee
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 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
    
    # solve the equlibrium
    d = K\F;
    
    return d,K
end
    # =======================================================================
    
    
    # =======================================================================
    
function  shape_quad4(r,s,xn,ien,nen, Imx_nsd,zero_nsd)
    
    #  computes the shape functions for quad4 element e at gauss point (r,s)
    # -----------------------------------------------------------------------
    # r               = gauss point coordinate in parent domain
    # s               = gauss point coordinate in parent domain
    # xn(nsd,nnp)     = nodal coordinates 
    # ien(nen,1)      = element connectivity for element e
    # nen             = number of nodes per element
    # Imx             = eye(nsd) 
    # Jaco            = zeros(nsd)
    #
    # sN(nen,1)       = shape functions at current gauss point
    # dNdx(nen,1)     = shape function derivatives at current gauss point in
    #                   (x,y)
    # dNdy(nen,1)     = shape function derivatives at current gauss point in
    #                   (x,y)
    # detJ            = determinant of the Jacobian for current gauss point
    # InvJ            = inverse of the Jacobian matrix for current gauss point
    #
    # =======================================================================
    
    
    # Initialize
    dNdx  = zeros(nen,1);
    dNdy  = zeros(nen,1);
    dNdy  .= dNdx;
    
    # compute shape functions at(r,s)
    sN,dNdr,dNds = shapefunct_quad4(r,s);
    
    # compute Jacobian, its deterninant and inverse
    Jaco,detJ,invJ = jacobian_2d(dNdr,dNds,xn[:,ien],nen,Imx_nsd,zero_nsd);
    
    # compute the derivatives of the shape functions
    for i = 1:nen
        dNdx[i] = invJ[1,1]*dNdr[i]+invJ[1,2]*dNds[i];
        dNdy[i] = invJ[2,1]*dNdr[i]+invJ[2,2]*dNds[i];
    end

    
    return sN,dNdx,dNdy,detJ,invJ
end
    # =======================================================================
    
    
    # =======================================================================
    
function  shapefunct_quad4(r,s)
    
    #  evaluates the shape functions and their derivatives at point (s,r)
    # -----------------------------------------------------------------------
    # r               = gauss point coordinate in parent domain
    # s               = gauss point coordinate in parent domain
    #
    # sN(nen,1)       = shape functions at current gauss point
    # dNdr(nen,1)     = shape function derivatives at current gauss point in
    #                   (r,s)
    # dNds(nen,1)     = shape function derivatives at current gauss point in
    #                   (r,s)
    #
    # =======================================================================
    
    # Initialize
    sN   = zeros(4,1);
    dNdr = zeros(4,1);
    dNds = zeros(4,1);

    dNdr .= sN;
    dNds .= sN;
    
    # Shape functions
    sN[1] = (1-r)*(1-s);
    sN[2] = (1+r)*(1-s);
    sN[3] = (1+r)*(1+s);
    sN[4] = (1-r)*(1+s);
    sN    .= 0.25*sN;
    
    # Derivatives
    dNdr[1] = -(1-s);
    dNdr[2] =  (1-s);
    dNdr[3] =  (1+s);
    dNdr[4] = -(1+s);
    dNdr    .= 0.25*dNdr;
    
    dNds[1] = -(1-r);
    dNds[2] = -(1+r);
    dNds[3] =  (1+r);
    dNds[4] =  (1-r);
    dNds    .= 0.25*dNds;
    
    return sN,dNdr,dNds
end
    # =======================================================================
    
# =======================================================================
function post_processing(d,E,g,id,ien,Ke,ndf,nee,nel,nen,nnp)

    # function that performs post processing for truss 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(1*nen,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];
    end

    # display("Fe")
    # display(Fe)
    return de,Fe
end