stiffnessMatrix.asv

% Written By: Matthew Jon Pais, University of Florida (2010)
% Website: http://sites.google.com/site/matthewjpais/Home
% Email: mpais@ufl.edu, matthewjpais@gmail.com

function [globalK, globalF] = stiffnessMatrix(omega,globalDOF,iter,enrElem,IElem)
% This function calculates the global stiffness matrix for the desired 
% discontinuities defined by the user supplied input.

global CHI CONNEC CRACK DOMAIN MAT NODES PHI PSI XYZ ZETA FORCE

nXElem  = DOMAIN(1);                                                        % Number of elements in the x-direction
nYElem  = DOMAIN(2);                                                        % Number of elements in the y-direction
lXElem  = DOMAIN(3);                                                        % Length of elements in the x-direction
Em      = MAT(1);                                                           % Young's modulus for the matrix
vm      = MAT(2);                                                           % Poisson's ratio for the matrix
Ef      = MAT(3);                                                           % Young's modulus for the fiber
vf      = MAT(4);                                                           % Poisson's ratio for the fiber
plane   = MAT(5);                                                           % Plane stress or plane strain
nCT     = size(PHI,2);                                                      % Number of crack tips 

ngausspt = 6; %number of Gauss points for quad integration cells
ngausspttri = 6; %number of Gauss points for tri integration cells
ngausspttip = 7; %number of Gauss points for tip integration cells

% Initialize the FE stiffness matrix
if iter > 1, globalDOF = globalDOF+16*iter; end                             % Initialize extra space for growing
globalK = sparse(globalDOF,globalDOF);                                      % Define the global K
globalF = zeros(1,globalDOF);
% Create elastic constant matrix
if plane == 1                                                               % Plane stress
    h = MAT(6);                                                             % Plane stress thickness
    
    C1 = Em/(1-vm^2);                                                       % Constant for elastic constant matrix
    C2 = Em*vm/(1-vm^2);                                                    % Constant for elastic constant matrix
    C3 = Em/2/(1+vm);                                                       % Constant for elastic constant matrix
    Cm = h*[C1 C2  0;...
            C2 C1  0;...
             0  0 C3];
    C1 = Ef/(1-vf^2);                                                       % Constant for elastic constant matrix
    C2 = Ef*vf/(1-vf^2);                                                    % Constant for elastic constant matrix
    C3 = Ef/2/(1+vf);                                                       % Constant for elastic constant matrix
    Cf = h*[C1 C2  0;...
            C2 C1  0;...
             0  0 C3];        
        
    % Constants for stiffness of matrix
    c0 = h*Em/(1-vm^2);
    c1 = c0*(1/2-vm/6);
    c2 = c0*(1/8+vm/8);
    c3 = c0*(-1/4-vm/12);
    c4 = c0*(-1/8+3*vm/8);
    c5 = c0*(-1/4+vm/12);
    c6 = c0*(-1/8-vm/8);
    c7 = c0*(vm/6);
    c8 = c0*(1/8-3*vm/8);

    % Stiffness of matrix
    Km = [c1 c2 c3 c4 c5 c6 c7 c8;
          c2 c1 c8 c7 c6 c5 c4 c3;
          c3 c8 c1 c6 c7 c4 c5 c2;
          c4 c7 c6 c1 c8 c3 c2 c5;
          c5 c6 c7 c8 c1 c2 c3 c4;
          c6 c5 c4 c3 c2 c1 c8 c7;
          c7 c4 c5 c2 c3 c8 c1 c6;
          c8 c3 c2 c5 c4 c7 c6 c1];
    
    % Constants for stiffness of fiber
    c0 = h*Ef/(1-vf^2);
    c1 = c0*(1/2-vf/6);
    c2 = c0*(1/8+vf/8);
    c3 = c0*(-1/4-vf/12);
    c4 = c0*(-1/8+3*vf/8);
    c5 = c0*(-1/4+vf/12);
    c6 = c0*(-1/8-vf/8);
    c7 = c0*(vf/6);
    c8 = c0*(1/8-3*vf/8);

    % Stiffness of fiber
    Kf = [c1 c2 c3 c4 c5 c6 c7 c8;
          c2 c1 c8 c7 c6 c5 c4 c3;
          c3 c8 c1 c6 c7 c4 c5 c2;
          c4 c7 c6 c1 c8 c3 c2 c5;
          c5 c6 c7 c8 c1 c2 c3 c4;
          c6 c5 c4 c3 c2 c1 c8 c7;
          c7 c4 c5 c2 c3 c8 c1 c6;
          c8 c3 c2 c5 c4 c7 c6 c1];         
elseif plane == 2                                                           % Plane strain
    C1 = Em*(1-vm)/(1+vm)/(1-2*vm);                                         % Constant for elastic constant matrix
    C2 = Em*vm/(1+vm)/(1-2*vm);                                             % Constant for elastic constant matrix
    C3 = Em/2/(1+vm);                                                       % Constant for elastic constant matrix
    Cm  = [C1 C2  0;...
           C2 C1  0;...
            0  0 C3];
    C1 = Ef*(1-vf)/(1+vf)/(1-2*vf);                                         % Constant for elastic constant matrix
    C2 = Ef*vf/(1+vf)/(1-2*vf);                                             % Constant for elastic constant matrix
    C3 = Ef/2/(1+vf);                                                       % Constant for elastic constant matrix
    Cf  = [C1 C2  0;...
           C2 C1  0;...
            0  0 C3];    
    
    % Constants for stiffness of matrix
    c0 = Em/(1+vm)/(1-2*vm);
    c1 = c0*(1/2-2*vm/3);
    c2 = c0*(1/8);
    c3 = c0*(-1/4+vm/6);
    c4 = c0*(-1/8+vm/2);
    c5 = c0*(-1/4+vm/3);
    c6 = c0*(-1/8);
    c7 = c0*(vm/6);
    c8 = c0*(1/8-vm/2);

    % Stiffness of matrix
    Km = [c1 c2 c3 c4 c5 c6 c7 c8;
          c2 c1 c8 c7 c6 c5 c4 c3;
          c3 c8 c1 c6 c7 c4 c5 c2;
          c4 c7 c6 c1 c8 c3 c2 c5;
          c5 c6 c7 c8 c1 c2 c3 c4;
          c6 c5 c4 c3 c2 c1 c8 c7;
          c7 c4 c5 c2 c3 c8 c1 c6;
          c8 c3 c2 c5 c4 c7 c6 c1];
    
    % Constants for stiffness of fiber
    c0 = Ef/(1+vf)/(1-2*vf);
    c1 = c0*(1/2-2*vf/3);
    c2 = c0*(1/8);
    c3 = c0*(-1/4+vf/6);
    c4 = c0*(-1/8+vf/2);
    c5 = c0*(-1/4+vf/3);
    c6 = c0*(-1/8);
    c7 = c0*(vf/6);
    c8 = c0*(1/8-vf/2);

    % Stiffness of fiber
    Kf = [c1 c2 c3 c4 c5 c6 c7 c8;
          c2 c1 c8 c7 c6 c5 c4 c3;
          c3 c8 c1 c6 c7 c4 c5 c2;
          c4 c7 c6 c1 c8 c3 c2 c5;
          c5 c6 c7 c8 c1 c2 c3 c4;
          c6 c5 c4 c3 c2 c1 c8 c7;
          c7 c4 c5 c2 c3 c8 c1 c6;
          c8 c3 c2 c5 c4 c7 c6 c1];   
end

m = size(CRACK,1);                                                          % Determine number of data points defining crack
if m > 0
    if nCT == 1
        xCT = CRACK(m,1);                                                   % X-coordinate of crack tip
        yCT = CRACK(m,2);                                                   % Y-coordinate of crack tip
    elseif nCT == 2
        xCT = [CRACK(1,1) CRACK(m,1)];                                      % X-coordinates of crack tips
        yCT = [CRACK(1,2) CRACK(m,2)];                                      % Y-coordinates of crack tips
    end
end

% Initialize the variables which will create the traditional sparse matrix
nIndexT  = 0;                                                               % Initialize index
uenrElem = nXElem*nYElem-length(enrElem)-length(IElem);                     % Number of unenriched elements
allRowT  = ones(uenrElem*64,1);                                             % Row indices
allColT  = ones(uenrElem*64,1);                                             % Column indices
allValT  = zeros(uenrElem*64,1);                                            % Stiffness matrix values

for iElem = 1:(nYElem*nXElem)
    N1  = CONNEC(iElem,2);                                                  % Node 1 for current element
    N2  = CONNEC(iElem,3);                                                  % Node 2 for current element
    N3  = CONNEC(iElem,4);                                                  % Node 3 for current element
    N4  = CONNEC(iElem,5);                                                  % Node 4 for current element
    NN  = NODES([N1 N2 N3 N4]',:);                                          % Nodal data for current element
    CN  = [CHI(N1) CHI(N2) CHI(N3) CHI(N4)];                                % Nodal chi level set values            
    CTN = nnz(NN(:,4));                                                     % Number of nodes with crack tip enrichment    
    HEN = nnz(NN(:,2));                                                     % Number of nodes with Heaviside enrichment
    IEN = nnz(NN(:,30));                                                    % Number of inclusion nodes    
    NEN = HEN+CTN+IEN;                                                      % Number of enriched nodes
  
    localK = 0;                                                             % Initialize stiffness for current element
    local  = [N1*2-1 N1*2 N2*2-1 N2*2 N3*2-1 N3*2 N4*2-1 N4*2];             % Traditional index locations
    iLoc   = 9;                                                             % Next index location
    localF=zeros(1, 8+2*NEN);  %size of local vector should be 8 + 2*number of enrichments
    if (NEN == 0)                                                           % Unenriched nodes
        ZN = [ZETA(N1) ZETA(N2) ZETA(N3) ZETA(N4)];
        if min(CN) >= 0                                                     % Traditional element
            if max(ZN) >= 0
                localK = Km;
            elseif max(ZN) < 0
                localK = Kf;
            end
        %elseif (max(CN) > 0) && (min(CN) < 0)
                       
            
            X1 = XYZ(N1,2); X2 = XYZ(N2,2); X3 = XYZ(N3,2); X4 = XYZ(N4,2); % Nodal x-coordinates
            Y1 = XYZ(N1,3); Y2 = XYZ(N2,3); Y3 = XYZ(N3,3); Y4 = XYZ(N4,3); % Nodal y-coordinates
            
            xyz = [X1 Y1;X2 Y2;X3 Y3;X4 Y4];                                % Nodal coordinate matrix            
            ZN  = [ZETA(N1) ZETA(N2) ZETA(N3) ZETA(N4)];                    % Nodal inclusion level set values
            [gp,gw] = gauss(ngausspt,'QUAD');
            %[gp,gw,J] = subDomain(3,CN,[],[],xyz,0,0,[]);
            xi        = gp(:,1);
            eta       = gp(:,2);
            N         = 1/4*[(1-xi).*(1-eta) (1+xi).*(1-eta) (1+xi).*(1+eta) (1-xi).*(1+eta)];
            chi       = N(:,1)*CN(1)+N(:,2)*CN(2)+N(:,3)*CN(3)+N(:,4)*CN(4);
            R         = find(chi <= 0);
            %J(R,:)    = [];
            %gp(R,:)   = [];
            %gw(R,:)   = [];
            
            %iElem
            %[gp, gw]
            
            for i = 1:size(gp,1)
                
                xi = gp(i,1); eta = gp(i,2);                                % Gauss points
                W  = gw(i);                                                 % Gauss weights
                
                %Ji   = [J(i,1) J(i,2);J(i,3) J(i,4)];                       % Jacobian of subdomain
                %detJ = det(Ji);                                             % Determinant of the Jacobian
                
                detJ  = lXElem/2*lXElem/2;
                
              %  Nx = 2/lXElem*1/4*[-(1-eta);1-eta;1+eta;-(1+eta)];          % Derivative of shape functions with respect to x
              %  Ny = 2/lXElem*1/4*[-(1-xi);-(1+xi);1+xi;1-xi];              % Derivative of shape functions with respect to y
                
                Zgp = N(i,1)*ZN(1)+N(i,2)*ZN(2)+N(i,3)*ZN(3)+N(i,4)*ZN(4);  % Material level set at current gauss point
                Xgp = N(i,1)*X1+N(i,2)*X2+N(i,3)*X3+N(i,4)*X4;
                Ygp = N(i,1)*Y1+N(i,2)*Y2+N(i,3)*Y3+N(i,4)*Y4;
                
                [fx, fy] = getBodyForce(Xgp,Ygp);
                
                localF(1:2:end) = localF(1:2:end) + W*N(i,:)*fx*detJ;
                localF(2:2:end) = localF(2:2:end) + W*N(i,:)*fy*detJ;
                                
             %   Bu = [Nx(1)   0   Nx(2)   0   Nx(3)   0   Nx(4)   0;...
             %           0   Ny(1)   0   Ny(2)   0   Ny(3)   0   Ny(4);...
             %         Ny(1) Nx(1) Ny(2) Nx(2) Ny(3) Nx(3) Ny(4) Nx(4)];
                  
                if Zgp > 0, C = Cm; else C = Cf; end
               
%                 C
                 localK = localK + W*Bu'*C*Bu*detJ;
%                 W
%                 Bu'
%                 C
%                 detJ
%                 pause
            end
%             iElem
%             Zgp
%             localK
%             localF
             %pause
            %localF
        end        
    elseif NEN > 0                                                          % Enriched element
        X1 = XYZ(N1,2); X2 = XYZ(N2,2); X3 = XYZ(N3,2); X4 = XYZ(N4,2);     % Nodal x-coordinates
        Y1 = XYZ(N1,3); Y2 = XYZ(N2,3); Y3 = XYZ(N3,3); Y4 = XYZ(N4,3);     % Nodal y-coordinates

        if NEN == 4                                                         % Fully enriched element
            xyz = [X1 Y1;X2 Y2;X3 Y3;X4 Y4];                                % Nodal coordinate matrix
            
            if numel(PSI) == 0, PN = [0 0 0 0]; else
                PN = [ PSI(N1)  PSI(N2)  PSI(N3)  PSI(N4)];                 % Nodal crack level set values
            end
            
            ZN = [ZETA(N1) ZETA(N2) ZETA(N3) ZETA(N4)];                     % Nodal inclusion level set values
            
            if IEN == 4
                [gp,gw,J] = subDomain(ngausspttri,CN,PN,ZN,xyz,0,0,[]);    
               
                
            elseif HEN == 4                                                 % Full Heaviside enrichment
                [gp,gw,J] = subDomain(ngausspttri,CN,PN,ZN,xyz,0,0,[]);
            elseif CTN == 4                                                 % Full crack tip enrichmnet
                [gp,gw,J] = subDomain(ngausspttip,CN,PN,ZN,xyz,1,0,[]);
            else                                                            % Full heaviside/crack tip enrichment
                [gp,gw,J] = subDomain(ngausspttip,CN,PN,ZN,xyz,0,0,[]);
              
            end
          
            if min(CN) < 0
                xi      = gp(:,1); 
                eta     = gp(:,2);                
                N       = 1/4*[(1-xi).*(1-eta) (1+xi).*(1-eta) (1+xi).*(1+eta) (1-xi).*(1+eta)];
                chi     = N(:,1)*CN(1)+N(:,2)*CN(2)+N(:,3)*CN(3)+N(:,4)*CN(4);
                R       = find(chi <= 0);
                J(R,:)  = [];
                gp(R,:) = [];
                gw(R,:) = [];
            end
          
        else                                                                % Partially enriched element
            if min(CN) < 0
                PN        = [ PSI(N1)  PSI(N2)  PSI(N3)  PSI(N4)];          % Nodal crack level set values
                ZN        = [ZETA(N1) ZETA(N2) ZETA(N3) ZETA(N4)];          % Nodal inclusion level set values
                
                [gp,gw,J] = subDomain(ngausspttri,CN,PN,ZN,xyz,0,0,[]);
                xi        = gp(:,1); 
                eta       = gp(:,2);                
                N         = 1/4*[(1-xi).*(1-eta) (1+xi).*(1-eta) (1+xi).*(1+eta) (1-xi).*(1+eta)];
                chi       = N(:,1)*CN(1)+N(:,2)*CN(2)+N(:,3)*CN(3)+N(:,4)*CN(4);
                R         = find(chi <= 0);
                J(R,:)  = [];
                gp(R,:) = [];
                gw(R,:) = [];              
            else
                [gp,gw] = gauss(ngausspt,'QUAD');
                J = [];
            end            
        end
        
        %gp
        %gw
       % NN
        for i = 1:length(gp)
            xi = gp(i,1); eta = gp(i,2);                                    % Gauss points
            W  = gw(i);                                                     % Gauss weights
            if isempty(J) == 0
                Ji   = [J(i,1) J(i,2);J(i,3) J(i,4)];                       % Jacobian of subdomain
                detJ = det(Ji);                                             % Determinant of the Jacobian
                
            else
                detJ  = lXElem/2*lXElem/2;                                  % Determinant of the Jacobian
            end

            N  = 1/4*[(1-xi)*(1-eta);(1+xi)*(1-eta);...                     % Shape functions
                      (1+xi)*(1+eta);(1-xi)*(1+eta)];
            Nx = 2/lXElem*1/4*[-(1-eta);1-eta;1+eta;-(1+eta)];              % Derivative of shape functions with respect to x
            Ny = 2/lXElem*1/4*[-(1-xi);-(1+xi);1+xi;1-xi];                  % Derivative of shape functions with respect to y

            X1 = XYZ(N1,2); X2 = XYZ(N2,2); X3 = XYZ(N3,2); X4 = XYZ(N4,2); % Nodal x-coordinates
            Y1 = XYZ(N1,3); Y2 = XYZ(N2,3); Y3 = XYZ(N3,3); Y4 = XYZ(N4,3); % Nodal y-coordinates
 
            Xgp = N(1)*X1+N(2)*X2+N(3)*X3+N(4)*X4;                          % The global X for the current gauss point
            Ygp = N(1)*Y1+N(2)*Y2+N(3)*Y3+N(4)*Y4;                          % The global Y for the current gauss point
            Zgp = N(1)*ZETA(N1)+N(2)*ZETA(N2)+N(3)*ZETA(N3)+N(4)*ZETA(N4);  % Material level set at current gauss point
            
            
            Benr = [];
            Nenr = [];
            Bu = [Nx(1)   0   Nx(2)   0   Nx(3)   0   Nx(4)   0;...
                    0   Ny(1)   0   Ny(2)   0   Ny(3)   0   Ny(4);...
                  Ny(1) Nx(1) Ny(2) Nx(2) Ny(3) Nx(3) Ny(4) Nx(4)];

            index = 1;
            
            for iN = 1:4
                if NN(iN,2) ~= 0
                    psi1 = PSI(N1);                                         % Psi level set value at node 1
                    psi2 = PSI(N2);                                         % Psi level set value at node 2
                    psi3 = PSI(N3);                                         % Psi level set value at node 3
                    psi4 = PSI(N4);                                         % Psi level set value at node 4
                    psi  = N(1)*psi1+N(2)*psi2+N(3)*psi3+N(4)*psi4;         % Psi level set value at current gauss point
    
                    Hgp = sign(psi);                                        % Heaviside value at current gauss point
                    Hi  = NN(iN,3);                                         % Nodal Heaviside value
                    H   = Hgp-Hi;                                           % Shifted Heaviside value

                    Ba = [Nx(iN)*H     0;
                              0    Ny(iN)*H;
                          Ny(iN)*H Nx(iN)*H];
                    Benr(:,index:(index+1)) = Ba;
                    index = index+2;
                    
                    if (i == length(gp))
                        local(iLoc:(iLoc+1)) = [2*NN(iN,2)-1 2*NN(iN,2)];
                        iLoc = iLoc+2;
                    end
                elseif NN(iN,4) ~= 0
                    if nCT == 1
                        X     = Xgp-xCT;                                    % Horizontal distance from crack tip to gauss point
                        Y     = Ygp-yCT;                                    % Vertical distance from crack tip to gauss point
                        CCS   = [cos(omega) sin(omega);-sin(omega) cos(omega)];
                        XYloc = CCS*[X Y]';                                 % Change to crack tip coordinates
                        r     = sqrt(XYloc(1)^2+XYloc(2)^2);                % Radius from crack tip to current gauss point
                        if r < 0.001*lXElem; r = 0.05*lXElem; end
                        theta = atan2(XYloc(2),XYloc(1));                   % Angle from crack tip to current gauss point
                    elseif nCT == 2
                        X1  = Xgp-xCT(1);
                        Y1  = Ygp-yCT(1);
                        X2  = Xgp-xCT(2);
                        Y2  = Ygp-yCT(2);
                        CCS = [cos(omega(1)) sin(omega(1));-sin(omega(1)) cos(omega(1))];
                        XY1 = CCS*[X1 Y1]';
                        CCS = [cos(omega(2)) sin(omega(2));-sin(omega(2)) cos(omega(2))];
                        XY2 = CCS*[X2 Y2]';
                        r1  = sqrt(XY1(1)^2+XY1(2)^2);                      % Radius from crack tip to current gauss point
                        r2  = sqrt(XY2(1)^2+XY2(2)^2);
                        if r1 > r2
                            r = r2; theta = atan2(XY2(2),XY2(1));
                            CCS = [cos(omega(2)) sin(omega(2));-sin(omega(2)) cos(omega(2))];                            
                        elseif r2 > r1
                            r = r1; theta = atan2(XY1(2),XY1(1));
                            CCS = [cos(omega(1)) sin(omega(1));-sin(omega(1)) cos(omega(1))];
                        end
                        if r < 0.001*lXElem; r = 0.05*lXElem; end     
                    end

                    c = 1/2/sqrt(r); ct = CCS(1,1); st = CCS(1,2);          % Constants
                                       
                    if NN(iN,12) == 0                                       % Crack tip enrichment
                        a1gp = sqrt(r)*sin(theta/2);                        % Node 1 crack tip enrichment value
                        a2gp = sqrt(r)*cos(theta/2);                        % Node 2 crack tip enrichment value
                        a3gp = sqrt(r)*sin(theta)*sin(theta/2);             % Node 3 crack tip enrichment value
                        a4gp = sqrt(r)*sin(theta)*cos(theta/2);             % Node 4 crack tip enrichment value

                        a1 = a1gp-NN(iN,5);                                 % Shifted alpha 1 enrichment value
                        a2 = a2gp-NN(iN,7);                                 % Shifted alpha 2 enrichment value
                        a3 = a3gp-NN(iN,9);                                 % Shifted alpha 3 enrichment value
                        a4 = a4gp-NN(iN,11);                                % Shifted alpha 4 enrichment value

                        % Derivative of crack tip enrichment functions with respect to X
                        Px = c*[-sin(theta/2)*ct              + cos(theta/2)*-st;...
                                 cos(theta/2)*ct              + sin(theta/2)*-st;...
                                -sin(3*theta/2)*sin(theta)*ct + (sin(theta/2)+sin(3*theta/2)*cos(theta))*-st;...
                                -cos(3*theta/2)*sin(theta)*ct + (cos(theta/2)+cos(3*theta/2)*cos(theta))*-st];

                        % Derivative of crack tip enrichment functions with respect to Y
                        Py = c*[-sin(theta/2)*st              + cos(theta/2)*ct;...
                                 cos(theta/2)*st              + sin(theta/2)*ct;...
                                -sin(3*theta/2)*sin(theta)*st + (sin(theta/2)+sin(3*theta/2)*cos(theta))*ct;...
                                -cos(3*theta/2)*sin(theta)*st + (cos(theta/2)+cos(3*theta/2)*cos(theta))*ct];

                        B1 = [Nx(iN)*a1+N(iN)*Px(1)          0;
                                       0            Ny(iN)*a1+N(iN)*Py(1);
                              Ny(iN)*a1+N(iN)*Py(1) Nx(iN)*a1+N(iN)*Px(1)];

                        B2 = [Nx(iN)*a2+N(iN)*Px(2)          0;
                                       0            Ny(iN)*a2+N(iN)*Py(2);
                              Ny(iN)*a2+N(iN)*Py(2) Nx(iN)*a2+N(iN)*Px(2)];

                        B3 = [Nx(iN)*a3+N(iN)*Px(3)          0;
                                       0            Ny(iN)*a3+N(iN)*Py(3);
                              Ny(iN)*a3+N(iN)*Py(3) Nx(iN)*a3+N(iN)*Px(3)];

                        B4 = [Nx(iN)*a4+N(iN)*Px(4)          0;
                                       0            Ny(iN)*a4+N(iN)*Py(4);
                              Ny(iN)*a4+N(iN)*Py(4) Nx(iN)*a4+N(iN)*Px(4)];

                        Bb = [B1 B2 B3 B4];
                        Benr(:,index:(index+7)) = Bb;
                        index = index+8;

                        if (i == length(gp))
                            local(iLoc:(iLoc+7)) = [2*NN(iN,4)-1 2*NN(iN,4) 2*NN(iN,6)-1  2*NN(iN,6)...
                                2*NN(iN,8)-1 2*NN(iN,8) 2*NN(iN,10)-1 2*NN(iN,10)];
                            iLoc = iLoc+8;
                        end
                    else                                                    % Bimaterial crack tip enrichment
                        G1 = Em/2/(1+vm);                                   % Shear modulus for the matrix
                        G2 = Ef/2/(1+vf);                                   % Shear modulus for the fiber                        
                        
                        % Kosolov constant
                        if plane == 1                                       % Plane stress
                            k1 = (3-vm)/(1+vm);
                            k2 = (3-vf)/(1+vf);
                        elseif plane == 2                                   % Plane strain
                            k1 = 3-4*vm;
                            k2 = 3-4*vf;
                        end
                        
                        b = (G1*(k2-1)-G2*(k1-1))/(G1*(k2+1)+G2*(k1+1));    % Second Dundur's parameter
                        e = 1/(2*pi)*log((1-b)/(1+b));                      % Material constant 
                        
                        % Common variables
                        sr  = sqrt(r);
                        st  = sin(theta);
                        st2 = sin(theta/2);
                        ct2 = cos(theta/2);
                       
                        a1gp  = sr*cos(e*log(r))*exp(-e*theta)*st2;         % Alpha 1 crack tip enrichment value
                        a2gp  = sr*cos(e*log(r))*exp(-e*theta)*ct2;         % Alpha 2 crack tip enrichment value
                        a3gp  = sr*cos(e*log(r))*exp(e*theta)*st2;          % Alpha 3 crack tip enrichment value
                        a4gp  = sr*cos(e*log(r))*exp(e*theta)*ct2;          % Alpha 4 crack tip enrichment value
                        a5gp  = sr*cos(e*log(r))*exp(e*theta)*st2*st;       % Alpha 5 crack tip enrichment value
                        a6gp  = sr*cos(e*log(r))*exp(e*theta)*ct2*st;       % Alpha 6 crack tip enrichment value
                        a7gp  = sr*sin(e*log(r))*exp(-e*theta)*st2;         % Alpha 7 crack tip enrichment value
                        a8gp  = sr*sin(e*log(r))*exp(-e*theta)*ct2;         % Alpha 8 crack tip enrichment value
                        a9gp  = sr*sin(e*log(r))*exp(e*theta)*st2;          % Alpha 9 crack tip enrichment value
                        a10gp = sr*sin(e*log(r))*exp(e*theta)*ct2;          % Alpha 10 crack tip enrichment value
                        a11gp = sr*sin(e*log(r))*exp(e*theta)*st2*st;       % Alpha 11 crack tip enrichment value
                        a12gp = sr*sin(e*log(r))*exp(e*theta)*ct2*st;       % Alpha 12 crack tip enrichment value
                        
                        a = [a1gp-NN(iN,5);...                              % Shifted alpha 1 enrichment value
                             a2gp-NN(iN,7);...                              % Shifted alpha 2 enrichment value
                             a3gp-NN(iN,9);...                              % Shifted alpha 3 enrichment value
                             a4gp-NN(iN,11);...                             % Shifted alpha 4 enrichment value
                             a5gp-NN(iN,13);...                             % Shifted alpha 5 enrichment value
                             a6gp-NN(iN,15);...                             % Shifted alpha 6 enrichment value
                             a7gp-NN(iN,17);...                             % Shifted alpha 7 enrichment value
                             a8gp-NN(iN,19);...                             % Shifted alpha 8 enrichment value
                             a9gp-NN(iN,21);...                             % Shifted alpha 9 enrichment value
                             a10gp-NN(iN,23);...                            % Shifted alpha 10 enrichment value
                             a11gp-NN(iN,25);...                            % Shifted alpha 11 enrichment value
                             a12gp-NN(iN,27)];                              % Shifted alpha 12 enrichment value
                        
                        % Derivative of bimaterial crack tip enrichment functions with respect to x1 (crack tip coordinate system)
                        px = c*[-exp(-e*theta)*sin(theta/2)*(cos(e*log(r))+2*e*sin(e*log(r)-theta));...
                                 exp(-e*theta)*cos(theta/2)*(cos(e*log(r))-2*e*sin(e*log(r)-theta));...
                                -exp(e*theta)*sin(theta/2)*(cos(e*log(r))+2*e*sin(e*log(r)+theta));...
                                 exp(e*theta)*cos(theta/2)*(cos(e*log(r))-2*e*sin(e*log(r)+theta));...
                                -exp(e*theta)*sin(theta)*(cos(e*log(r))*sin(3*theta/2)+2*e*sin(e*log(r)+theta)*sin(theta/2));...
                                -exp(e*theta)*sin(theta)*(cos(e*log(r))*cos(3*theta/2)+2*e*sin(e*log(r)+theta)*cos(theta/2));...
                                 exp(-e*theta)*sin(theta/2)*(-sin(e*log(r))+2*e*cos(e*log(r)-theta));...
                                 exp(-e*theta)*cos(theta/2)*(sin(e*log(r))+2*e*cos(e*log(r)-theta));...
                                 exp(e*theta)*sin(theta/2)*(-sin(e*log(r))+2*e*cos(e*log(r)+theta));...
                                 exp(e*theta)*cos(theta/2)*(sin(e*log(r))+2*e*cos(e*log(r)+theta));...
                                 exp(e*theta)*sin(theta)*(-sin(e*log(r))*sin(3*theta/2)+2*e*cos(e*log(r)+theta)*sin(theta/2));...
                                 exp(e*theta)*sin(theta)*(-sin(e*log(r))*cos(3*theta/2)+2*e*cos(e*log(r)+theta)*cos(theta/2))];

                        % Derivative of bimaterial crack tip enrichment functions with respect to x2 (crack tip coordinate system)
                        py = c*[exp(-e*theta)*(cos(e*log(r))*cos(theta/2)-2*e*cos(e*log(r)-theta)*sin(theta/2));...
                                exp(-e*theta)*(cos(e*log(r))*sin(theta/2)-2*e*cos(e*log(r)-theta)*cos(theta/2));...
                                exp(e*theta)*(cos(e*log(r))*cos(theta/2)+2*e*cos(e*log(r)+theta)*sin(theta/2));...
                                exp(e*theta)*(cos(e*log(r))*sin(theta/2)+2*e*cos(e*log(r)+theta)*cos(theta/2));...
                                exp(e*theta)*(cos(e*log(r))*(sin(theta/2)+sin(3*theta/2)*cos(theta))+2*e*cos(e*log(r)+theta)*sin(theta/2)*sin(theta));...
                                exp(e*theta)*(cos(e*log(r))*(cos(theta/2)+cos(3*theta/2)*cos(theta))+2*e*cos(e*log(r)+theta)*cos(theta/2)*sin(theta));...
                                exp(-e*theta)*(sin(e*log(r))*cos(theta/2)-2*e*sin(e*log(r)-theta)*sin(theta/2));...
                                exp(-e*theta)*(sin(e*log(r))*sin(theta/2)-2*e*sin(e*log(r)-theta)*cos(theta/2));...
                                exp(e*theta)*(sin(e*log(r))*cos(theta/2)+2*e*sin(e*log(r)+theta)*sin(theta/2));...
                                exp(e*theta)*(sin(e*log(r))*sin(theta/2)+2*e*sin(e*log(r)+theta)*cos(theta/2));...
                                exp(e*theta)*(sin(e*log(r))*(sin(theta/2)+sin(3*theta/2)*cos(theta))+2*e*sin(e*log(r)+theta)*sin(theta/2)*sin(theta));...
                                exp(e*theta)*(sin(e*log(r))*(cos(theta/2)+cos(3*theta/2)*cos(theta))+2*e*sin(e*log(r)+theta)*cos(theta/2)*sin(theta))];
                        
                        % Derivative of bimaterial crack tip enrichment functions with respect to X
                        Px = [px(1)*ct+py(1)*-st;...
                              px(2)*ct+py(2)*-st;...
                              px(3)*ct+py(3)*-st;...
                              px(4)*ct+py(4)*-st;...
                              px(5)*ct+py(5)*-st;...
                              px(6)*ct+py(6)*-st;...
                              px(7)*ct+py(7)*-st;...
                              px(8)*ct+py(8)*-st;...
                              px(9)*ct+py(9)*-st;...
                              px(10)*ct+py(10)*-st;...
                              px(11)*ct+py(11)*-st;...
                              px(12)*ct+py(12)*-st];

                        % Derivative of bimaterial crack tip enrichment functions with respect to Y
                        Py = [px(1)*st+py(1)*ct;...
                              px(2)*st+py(2)*ct;...
                              px(3)*st+py(3)*ct;...
                              px(4)*st+py(4)*ct;...
                              px(5)*st+py(5)*ct;...
                              px(6)*st+py(6)*ct;...
                              px(7)*st+py(7)*ct;...
                              px(8)*st+py(8)*ct;...
                              px(9)*st+py(9)*ct;...
                              px(10)*st+py(10)*ct;...
                              px(11)*st+py(11)*ct;...
                              px(12)*st+py(12)*ct];

                        Bb = zeros(3,24);
                        for iB = 1:12
                            Balpha = [Nx(iN)*a(iB)+N(iN)*Px(iB)             0;...
                                                  0             Nx(iN)*a(iB)+N(iN)*Py(iB);...
                                      Nx(iN)*a(iB)+N(iN)*Py(iB) Nx(iN)*a(iB)+N(iN)*Px(iB)];
                            Bb(:,(2*iB-1):2*iB) = Balpha;
                        end
                        Benr(:,index:(index+23)) = Bb;
                        index = index+24;
                        if i == length(gp)
                            for iB = 4:2:26
                                local(iLoc:(iLoc+1)) = [2*NN(iN,iB)-1 2*NN(iN,iB)];
                                iLoc = iLoc+2;
                            end
                        end        
                    end
                end
                
                % Inclusion enrichment
                if NN(iN,30) ~= 0
                    zeta = ZETA([N1 N2 N3 N4]);
                    zabs = abs(zeta);
                    Zm   = N(1)*zeta(1)+N(2)*zeta(2)+N(3)*zeta(3)+N(4)*zeta(4);%dot(N,zeta);
                    Za   = N(1)*zabs(1)+N(2)*zabs(2)+N(3)*zabs(3)+N(4)*zabs(4);%dot(N,zabs);
                    Zmn  = Zm/abs(Zm);

                    %E  = Za-abs(Zm);
                    %E = abs(Zm)-zabs(iN); %Sukumar
                    %Ex = dot(Nx(1:4),zeta(1:4))*sign(Zm);
                    %Ey = dot(Ny(1:4),zeta(1:4))*sign(Zm);
                    

                    E = Za-abs(Zm); %Moes
                    %Ex = dot(Nx(1:4),zabs(1:4)) - dot(Nx(1:4),zeta(1:4))*sign(Zm);
                    %Ey = dot(Ny(1:4),zabs(1:4)) - dot(Ny(1:4),zeta(1:4))*sign(Zm);
                    
              %      E = 1;
              %      Ex = 0;
              %      Ey = 0;
                    
                     Ex = Nx(1)*zabs(1)+Nx(2)*zabs(2)+Nx(3)*zabs(3)+Nx(4)*zabs(4)-...
                          Zmn*(Nx(1)*zeta(1)+Nx(2)*zeta(2)+Nx(3)*zeta(3)+Nx(4)*zeta(4));%dot(Nx,zabs)-Zm/zmab*dot(Nx,zeta);
                     Ey = Ny(1)*zabs(1)+Ny(2)*zabs(2)+Ny(3)*zabs(3)+Ny(4)*zabs(4)-...
                          Zmn*(Ny(1)*zeta(1)+Ny(2)*zeta(2)+Ny(3)*zeta(3)+Ny(4)*zeta(4));%dot(Ny,zabs)-Zm/zmab*dot(Ny,zeta);
                                       
                    Ba = [Nx(iN)*E+N(iN)*Ex         0;
                                  0         Ny(iN)*E+N(iN)*Ey;
                          Ny(iN)*E+N(iN)*Ey Nx(iN)*E+N(iN)*Ex];
                    Benr(:,index:(index+1)) = Ba;
                    %iN
                    index = index+2;
                    %Nenr = [Nenr, N(iN)*E];
                    Nenr((index-1)/2) = N(iN)*E;
                    %Nenr((index-1)/2) = E;
                    %Nenr(index) = N(iN)*E;
                    %Nenr(index+1) = N(iN)*E;
                    if (i == length(gp))
                        local(iLoc:(iLoc+1)) = [2*NN(iN,30)-1 2*NN(iN,30)];
                        iLoc = iLoc+2;
                    end      
                    %length(gp)

                end
            end

            
            if Zgp > 0, C = Cm; else C = Cf; end
            B = [Bu Benr];            
            localK = localK + W*B'*C*B*detJ; 
            
            %N
            %Nenr
            
            %Xgp
            %Ygp
            Nenr = Nenr;
            Nall = [N' Nenr];
            
            %detJ
           
            %option 1: fx = FORCE(4); fy = FORCE(5);
          
            [fx, fy] = getBodyForce(Xgp,Ygp);
              %temp = W*Nall*fx*detJ;
              %endindex = length(temp);
              %W
              %Nall
              %sum(Nall(1:4))
              %fy
              %detJ
              %temp = W*Nall*fy*detJ;
              
              %pause
              
             
             localF(1:2:end) = localF(1:2:end) + W*Nall*fx*detJ;
             localF(2:2:end) = localF(2:2:end) + W*Nall*fy*detJ;
             %disp(detJ)
            % oddlocalF = localF(1:2:end-1)
            % evanlocalF = localF(2:2:end)
             
             %localF(1:2:end-1) = localF(1:2:end-1) + W*Nall(i,:)*fx*detJ;
             %localF(2:2:end) = localF(2:2:end) + W*Nall(i,:)*fy*detJ;
             %localF = localF + W*[Nall*fx, Nall*fy]'*detJ;
             %pause
             
             %disp(localF);
             %localF(2:2:end)
             
        end
        %pause
    end 
    
    if length(localK) == 8                                                  % Unenriched element
        
        for row = 1:8
            for col = 1:8
                nIndexT = nIndexT+1;
                allRowT(nIndexT) = local(row);
                allColT(nIndexT) = local(col);
                allValT(nIndexT) = localK(row,col);
            end
        end
        
        
    else
        %disp(local);
        %pause;
        %globalK(local,local) = globalK(local,local) + localK;               % Assemble the global stiffness
        %globalF(local) = globalF(local) + localF; 
        %disp(globalF)
    end
    %localF
    %localK
   % pause
    %disp(localF)
    globalK(local,local) = globalK(local,local) + localK;               % Assemble the global stiffness
    globalF(local) = globalF(local) + localF; 
    %size(globalF(local))
    %size(localF)
    %localF
    %pause
    %globalF(local) = globalF(local) + localF;
    %localF
    %pause
end

%globalK = globalK + sparse(allRowT,allColT,allValT,globalDOF,globalDOF);