center_crack_1.m

efg code with matlab

        for m = 1 : le
            B_fem = [dphidx(m) 0 ;
                     0 dphidy(m) ;
                     dphidy(m) dphidx(m)];
            % a split enriched node
            if (snode(m) ~= 0)
                % compute Heaviside and derivatives
                dist = signed_distance(xCr,pt);
                [H,dHdx,dHdy] = heaviside(dist);

                B_enr = [dphidx(m)*H 0 ;
                         0 dphidy(m)*H;
                         dphidy(m)*H dphidx(m)*H];
                Bm = [B_fem B_enr];     
                clear B_fem ; clear B_enr;
                % compute positions in K matrix
                en(2*ind+1) = u_pos(index(m))  ;
                en(2*ind+2) = u_pos(index(m))+1;
                en(2*ind+3) = u_pos(index(m))+2;
                en(2*ind+4) = u_pos(index(m))+3;
                ind = ind + 2 ;
            elseif (tnode(m) ~= 0)
                % compute branch functions
                xp = QT*(pt-xTip)';                % local coordinates
                [theta,r] = cart2pol(xp(1),xp(2)); % local polar coordinates
                [Br,dBdx,dBdy] = branch(r,theta,alpha);
                B1_enr = [dphidx(m)*Br(1) + phi(m)*dBdx(1) 0 ;
                          0 dphidy(m)*Br(1) + phi(m)*dBdy(1);
                          dphidy(m)*Br(1) + phi(m)*dBdy(1) dphidx(m)*Br(1) + phi(m)*dBdx(1)];

                B2_enr = [dphidx(m)*Br(2) + phi(m)*dBdx(2) 0 ;
                          0 dphidy(m)*Br(2) + phi(m)*dBdy(2);
                          dphidy(m)*Br(2) + phi(m)*dBdy(2) dphidx(m)*Br(2) + phi(m)*dBdx(2)];

                B3_enr = [dphidx(m)*Br(3) + phi(m)*dBdx(3) 0 ;
                          0 dphidy(m)*Br(3) + phi(m)*dBdy(3);
                          dphidy(m)*Br(3) + phi(m)*dBdy(3) dphidx(m)*Br(3) + phi(m)*dBdx(3)];

                B4_enr = [dphidx(m)*Br(4) + phi(m)*dBdx(4) 0 ;
                          0 dphidy(m)*Br(4) + phi(m)*dBdy(4);
                          dphidy(m)*Br(4) + phi(m)*dBdy(4) dphidx(m)*Br(4) + phi(m)*dBdx(4)];
                B_enr = [B1_enr B2_enr B3_enr B4_enr];
                clear B1_enr; clear B2_enr; clear B3_enr; clear B4_enr;
                Bm = [B_fem B_enr]; 
                clear B_enr; 
                % compute positions in K matrix
                en(2*ind+1)  = u_pos(index(m))  ;
                en(2*ind+2)  = u_pos(index(m))+1;
                en(2*ind+3)  = u_pos(index(m))+2;
                en(2*ind+4)  = u_pos(index(m))+3;
                en(2*ind+5)  = u_pos(index(m))+4;
                en(2*ind+6)  = u_pos(index(m))+5;
                en(2*ind+7)  = u_pos(index(m))+6;
                en(2*ind+8)  = u_pos(index(m))+7;
                en(2*ind+9)  = u_pos(index(m))+8;
                en(2*ind+10) = u_pos(index(m))+9;
                ind = ind + 5 ;
            else
                en(2*ind+1) = u_pos(index(m))   ;
                en(2*ind+2) = u_pos(index(m))+1 ;
                ind = ind + 1;
                Bm = B_fem;
            end  
            B = [B Bm];  
            clear B_fem;
        end   % end of loop on nodes in neighbour of Gp
        K(en,en) = K(en,en) + B'*C*B * W(igp)*J(igp) ;
    end       % end of check enriched nodes
end          % end of loop on Gauss Points

% +++++++++++++++++++++++++++++++++++++
%    INTEGRATE ON TRACTION BOUNDARY
% +++++++++++++++++++++++++++++++++++++
% Steps:
% 1. Build up Gauss points on the boundary, 4 GPs are used
% 2. Loop on these GPs to form the f vector
disp([num2str(toc),'   No imposed traction in this problem'])


% +++++++++++++++++++++++++++++++++++++
%   INTEGRATE ON DISPLACEMENT BOUNDARY
% +++++++++++++++++++++++++++++++++++++
% Steps:
% 1. Build up Gauss points on the boundary, 4 GPs are used
% 2. Loop on these GPs to form the vector q and matrix G

disp([num2str(toc),'   INTEGRATION ON DISPLACEMENT BOUNDARY'])

% If Lagrange multiplier method is used
if disp_bc_method == 1
    disp(' Lagrange multiplier method is used ')
    % --------------------------------------
    %       Form vector q and matrix G
    % --------------------------------------
    qk = zeros(1,2*num_disp_nodes);
    G  = zeros(2*total_num_node,2*num_disp_nodes);

    % If point collocation method is used. Assume that these collocation
    % points are coincident to nodes on displacement boundary
    if (la_approx == 200)
        disp(' Point collocation method is used ')
        for i = 1 : num_disp_nodes      % loop on collocation points
            pt = node(disp_nodes(i),:);
            % compute exact displacement, u_bar
            [ux,uy] = exact_Griffith(pt,E,nu,stressState,sigmato,xTip,seg,cracklength);
            
            % qk vector 
            qk(1,2*i-1) = -ux ;
            qk(1,2*i)   = -uy ;
            
            % S matrix
            S = [1 0; 0 1]; % y direction
            % G matrix
            [index] = define_support(node,pt,di);
            [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);
            for j = 1 : size(index,2)
                row1 = u_pos(index(j))     ; 
                row2 = u_pos(index(j)) + 1 ;
                Gj   =  -phi(j) * S ;
                G(row1:row2,2*i-1:2*i) = G(row1:row2,2*i-1:2*i) + Gj ;

            end
        end
    end   % end of if boundary point collocation is used

    % if Lagrange multipliers are interpolated by finite elements
    % along the displacement boundary
    if (la_approx == 100)
        disp(' Finite element interpolation for ML is used ')
        m1 = 0 ;
        for i = 1 : (num_disp_nodes - 1)
            sctr = [disp_nodes(i) disp_nodes(i+1)];
            m1 = m1 + 1 ;
            m2 = m1 + 1 ;
            le = norm(disp_nodes(m1) - disp_nodes(m2));
            for q = 1:size(W1,1)
                pt = Q1(q,:);
                wt = W1(q);
                [N,dNdxi] = lagrange_basis('L2',pt);
                J0 = dNdxi'*node(sctr,:);
                detJ = norm(J0) ;
                pt = N' * node(sctr,:); % global GP
                % compute exact displacement of Timoshenko problem, u_bar
                x = pt(1,1) ; y = pt(1,2) ;
                fac = P/(6*E*I);
                ux_exact =  fac*y*((6*L-3*x)*x+(2+nu)*(y^2-D*D/4));
                uy_exact = -fac*((L-x)*3*nu*y^2+0.25*(4+5*nu)*D*D*x+(3*L-x)*x^2);

                N1 = 1 - pt(1,2)/le ; N2 = 1 - N1 ;
                % qk vector
                qk(2*m1-1) = qk(2*m1-1) - wt * detJ * N1 * ux_exact ;
                qk(2*m1)   = qk(2*m1)   - wt * detJ * N1 * uy_exact ;
                qk(2*m2-1) = qk(2*m2-1) - wt * detJ * N2 * ux_exact ;
                qk(2*m2)   = qk(2*m2)   - wt * detJ * N2 * uy_exact ;

                % G matrix
                [index] = define_support(node,pt,di);
                [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);
                for j = 1 : size(index,2)
                    row1 = 2*index(j)-1 ;
                    row2 = 2*index(j)   ;
                    G1   = - wt * detJ * phi(j) * [N1 0 ; 0 N1];
                    G2   = - wt * detJ * phi(j) * [N2 0 ; 0 N2];
                    G(row1:row2,2*m1-1:2*m1) = G(row1:row2,2*m1-1:2*m1) + G1;
                    G(row1:row2,2*m2-1:2*m2) = G(row1:row2,2*m2-1:2*m2) + G2;
                end
            end
        end
    end
end       % end of if disp_bc_method = 'Lagrange'

% If penalty method is use, then one must modify the stiffness matrix and
% the nodal force vector
% $K_{ij}$ = Kij - alpha \int phi_i phi_j d \gamma_u
% fj  = fj - alpha \int phi_i u_bar d \gamma_u

if disp_bc_method == 2
    fu = zeros(2*numnode,1);
    phi_ij = 0;
    for igp = 1 : size(Wu,1)
        pt = Qu(igp,:);                             % quadrature point
        wt = Wu(igp);                               % quadrature weight
        [index] = define_support(node,pt,di);
        [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);
        % compute exact displacement of Timoshenko problem
        x = pt(1,1) ; y = pt(1,2) ;
        fac = P/(6*E*I);
        ux_exact =  fac*y*((6*L-3*x)*x+(2+nu)*(y^2-D*D/4));
        uy_exact = -fac*((L-x)*3*nu*y^2+0.25*(4+5*nu)*D*D*x+(3*L-x)*x^2);
        le = size(index,2);
        en = zeros(1,2*le);
        force = zeros(1,2*le);
        for j = 1 : le
            en(2*j-1) = 2*index(j)-1;
            en(2*j  ) = 2*index(j)  ;
            force(2*j-1) = phi(j)*ux_exact;
            force(2*j  ) = phi(j)*uy_exact;
        end
        fu(en) = fu(en) + Jt(igp) * wt * force' ;
        phi_ij = phi_ij + Jt(igp)*wt*phi*phi' ;
    end % end of loop on Gauss points
    f = f - alpha*fu;
    K(en,en) = K(en,en) - alpha*phi_ij ;
end

% +++++++++++++++++++++++++++++++++++++
%      SOLUTION OF THE EQUATIONS
% +++++++++++++++++++++++++++++++++++++

disp([num2str(toc),'   SOLUTION OF EQUATIONS'])

% If penalty function is used, then the unknowns u is the same, just solve
% to get it:-)
if disp_bc_method == 2
    d = K\f ;
    for i = 1 : numnode
        u(1,i) = d(2*i-1); % x displacement
        u(2,i) = d(2*i);   % y displacement
    end
end

% If Lagrange multiplier is used then, extend the unknowns u by lamda
if disp_bc_method == 1
    f = [f;qk'];                                % f = {f;qk}
    m = ([K G; G' zeros(2*num_disp_nodes)]);    % m = [K GG;GG' 0]
    d = m\f;                                    % d = {u;lamda}
    % just get displacement parameters u_i
    u = d(1:2*total_num_node);
end

clear d ;

% +++++++++++++++++++++++++++++++++++++
%    COMPUTE THE TRUE DISPLACEMENTS
% +++++++++++++++++++++++++++++++++++++

disp([num2str(toc),'   INTERPOLATION TO GET TRUE DISPLACEMENT'])

% +++++++++++++++++++++++++++++++++++++
%       COMPUTE STRESSES AT NODES
% +++++++++++++++++++++++++++++++++++++
disp([num2str(toc),'   COMPUTE STRESS AT GAUSS POINTS'])
count = 0 ;
for i = 1 : numnode
    pt = node(i,:);                              
    wt = W(igp);                              
    [index] = define_support(node,pt,di);
    le      = size(index,2);
    [snode,s_loc] = ismember(index,split_nodes);
    [tnode,t_loc] = ismember(index,tip_nodes);
    if (any(snode) == 0)
        no_split_node  = 1 ;
        num_split_node = 0 ;
    else
        no_split_node = 0;
        num_split_node = size(find(snode ~= 0),2);
    end
    if (any(tnode) == 0)
        no_tip_node = 1;
        num_tip_node = 0;
    else
        no_tip_node = 0;
        num_tip_node = size(find(tnode ~= 0),2);
    end
    if (no_tip_node == 1) && (no_split_node == 1)
        B  = zeros(3,2*le) ;
        en = zeros(1,2*le);
        [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);
        for m = 1 : le
            B(1:3,2*m-1:2*m) = [dphidx(m) 0 ;
                                0 dphidy(m);
                                dphidy(m) dphidx(m)];
            en(2*m-1) = u_pos(index(m))    ;
            en(2*m  ) = u_pos(index(m))+1  ;
        end
    else  % there are enriched nodes, attention !!!
        le_of_B = (le + num_split_node*1 + num_tip_node*4) ;
        en = zeros(1,2*le_of_B);
        B = []; 
        [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);
        ind = 0 ;
        for m = 1 : le
            B_fem = [dphidx(m) 0 ;
                     0 dphidy(m) ;
                     dphidy(m) dphidx(m)];
            % a split enriched node
            if (snode(m) ~= 0)
                disp('ah, one split enriched node')
                % compute Heaviside and derivatives
                dist = signed_distance(xCr,pt);
                [H,dHdx,dHdy] = heaviside(dist);

                B_enr = [dphidx(m)*H 0 ;
                         0 dphidy(m)*H;
                         dphidy(m)*H dphidx(m)*H];
                Bm = [B_fem B_enr];     
                clear B_fem ; clear B_enr;
                % compute positions in K matrix
                en(2*ind+1) = u_pos(index(m))  ;
                en(2*ind+2) = u_pos(index(m))+1;
                en(2*ind+3) = u_pos(index(m))+2;
                en(2*ind+4) = u_pos(index(m))+3;
                ind = ind + 2 ;
            elseif (tnode(m) ~= 0)
                disp('ah, one tip enriched node')
                % compute branch functions
                xp = QT*(pt-xTip)';                % local coordinates
                [theta,r] = cart2pol(xp(1),xp(2)); % local polar coordinates
                [Br,dBdx,dBdy] = branch(r,theta,alpha);
                B1_enr = [dphidx(m)*Br(1) + phi(m)*dBdx(1) 0 ;
                          0 dphidy(m)*Br(1) + phi(m)*dBdy(1);
                          dphidy(m)*Br(1) + phi(m)*dBdy(1) dphidx(m)*Br(1) + phi(m)*dBdx(1)];

                B2_enr = [dphidx(m)*Br(2) + phi(m)*dBdx(2) 0 ;
                          0 dphidy(m)*Br(2) + phi(m)*dBdy(2);
                          dphidy(m)*Br(2) + phi(m)*dBdy(2) dphidx(m)*Br(2) + phi(m)*dBdx(2)];

                B3_enr = [dphidx(m)*Br(3) + phi(m)*dBdx(3) 0 ;
                          0 dphidy(m)*Br(3) + phi(m)*dBdy(3);
                          dphidy(m)*Br(3) + phi(m)*dBdy(3) dphidx(m)*Br(3) + phi(m)*dBdx(3)];

                B4_enr = [dphidx(m)*Br(4) + phi(m)*dBdx(4) 0 ;
                          0 dphidy(m)*Br(4) + phi(m)*dBdy(4);
                          dphidy(m)*Br(4) + phi(m)*dBdy(4) dphidx(m)*Br(4) + phi(m)*dBdx(4)];
                B_enr = [B1_enr B2_enr B3_enr B4_enr];
                clear B1_enr; clear B2_enr; clear B3_enr; clear B4_enr;
                Bm = [B_fem B_enr]; 
                clear B_enr; 
                % compute positions in K matrix
                en(2*ind+1)  = u_pos(index(m))  ;
                en(2*ind+2)  = u_pos(index(m))+1;
                en(2*ind+3)  = u_pos(index(m))+2;
                en(2*ind+4)  = u_pos(index(m))+3;
                en(2*ind+5)  = u_pos(index(m))+4;
                en(2*ind+6)  = u_pos(index(m))+5;
                en(2*ind+7)  = u_pos(index(m))+6;
                en(2*ind+8)  = u_pos(index(m))+7;
                en(2*ind+9)  = u_pos(index(m))+8;
                en(2*ind+10) = u_pos(index(m))+9;
                ind = ind + 5 ;
            else
                en(2*ind+1) = u_pos(index(m))   ;
                en(2*ind+2) = u_pos(index(m))+1 ;
                ind = ind + 1;
                Bm = B_fem;
            end  
            B = [B Bm];  
            clear B_fem;
        end   
    end 
    count = count + 1 ;
    stress(1:3,count) = C*B*u(en);  % sigma = C*epsilon = C*(B*u)
end          

figure
plot_field(node,element,'Q4',stress(1,:));
axis('equal');
xlabel('X');
ylabel('Y');
title('Sigma XX');
set(gcf,'color','white');
colorbar('vert');