center_crack.m

efg code with matlab

% +++++++++++++++++++++++++++++++++++++
% 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 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]; % both directions
            % 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)    ;
                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 ')
        
        % 4 point quadrature for each "1D element"
        [W1,Q1]=quadrature(4, 'GAUSS', 1 );

        % ------------------------------
        %   Generation of Gauss points
        % ------------------------------
        Wu = [];
        Qu = [];
        Ju = [];

        for i = 1 : (num_disp_nodes - 1)
            sctr = [disp_nodes(i) disp_nodes(i+1)];
            for q = 1:size(W1,1)
                pt = Q1(q,:);
                wt = W1(q);
                [N,dNdxi] = lagrange_basis('L2',pt);
                J0 = dNdxi'*node(sctr,:);
                Qu = [Qu; N' * node(sctr,:)];
                Wu = [Wu; wt];
                Ju = [Ju; norm(J0)];
            end % of quadrature loop
        end
        
        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
                [ux,uy] = exact_Griffith(pt,E,nu,stressState,sigmato,xTip,seg,...
                    cracklength);
                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'


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

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


% 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(num_disp_nodes*2)]);    % m = [K GG;GG' 0]
    d = m\f;                                    % d = {u;lamda}
    % just get nodal parameters u_i, not need Lagrange multipliers
    u = d(1:2*total_num_node);
end

clear d ;

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

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

for i = 1 : numnode
    pt = node(i,:);
    [index] = define_support(node,pt,di);
    le      = length(index);

    [snode,s_loc] = ismember(index,split_nodes);
    [tnode,t_loc] = ismember(index,tip_nodes);

    if (any(snode) == 0)
        no_split_node  = 1 ; % no H(x) enriched node
        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;    % no tip enriched node
        num_tip_node = 0;
    else
        no_tip_node = 0;
        num_tip_node = size(find(tnode ~= 0),2);
    end

    % shape function at nodes in neighbouring of node i
    [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);

    % Recall the enriched displacement approximation:
    % u(x) = Ni*ui + Nj*H*bj + Nk*Bk*bk

    % The first part Ni*ui
    stdUx = phi*u(index.*2-1);
    stdUy = phi*u(index.*2);

    if (no_tip_node == 1) && (no_split_node == 1)
        ux(i) = stdUx ;
        uy(i) = stdUy ;
    else
        % The second and third part exist only with enriched particles
        HUx = 0 ;
        HUy = 0 ;
        BUx = 0 ;
        BUy = 0 ;
        for m = 1 : le
            % a split enriched node
            if (snode(m) ~= 0)
                % compute Heaviside and derivatives
                dist = signed_distance(xCr,pt);
                [H,dHdx,dHdy] = heaviside(dist);
                HUx = HUx + phi(m)*H*u(2 * pos(index(m)) - 1) ;
                HUy = HUy + phi(m)*H*u(2 * pos(index(m))    ) ;
            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);
                BUx = BUx + phi(m)*Br(1)*u(2 * pos(index(m)) - 1)+... 
                            phi(m)*Br(2)*u(2 * (pos(index(m))+1) - 1)+...
                            phi(m)*Br(3)*u(2 * (pos(index(m))+2) - 1)+...
                            phi(m)*Br(4)*u(2 * (pos(index(m))+3) - 1);         
                BUy = BUy + phi(m)*Br(1)*u(2 * pos(index(m)))+... 
                            phi(m)*Br(2)*u(2 * (pos(index(m))+1))+...
                            phi(m)*Br(3)*u(2 * (pos(index(m))+2))+...
                            phi(m)*Br(4)*u(2 * (pos(index(m))+3));                              
            end
        end   % end of loop on nodes in neighbour of Gp
        ux(i) = stdUx + HUx + BUx ;
        uy(i) = stdUy + HUy + BUy ;        
    end
end

% ++++++++++++++++++++++++++++++++++++
%    COMPUTE STRESSES AT PARTICLES
% ++++++++++++++++++++++++++++++++++++
disp([num2str(toc),'   STRESSES COMPUTATION'])

ind = 0 ;
for igp = 1 : size(node,1)
    pt = node(igp,:);
    [index] = define_support(node,pt,di);
    [sctrB,snode,tnode] = assembly(index,split_nodes,tip_nodes,...
                                        pos);
    %  compute B matrix
    B = Bmatrix(pt,index,node,di,form,...
                     snode,tnode,xCr,xTip,alpha);
    ind = ind + 1 ;
    stress_gp(1:3,ind) = C*B*u(sctrB);  % sigma = C*epsilon = C*(B*u)
end

fac = 30 ; % visualization factor
figure
plot_field(node+fac*[ux' uy'],element,'Q4',stress_gp(1,:));
axis('equal');
xlabel('X');
ylabel('Y');
title('\sigma_{xx}');
set(gcf,'color','white');
colorbar('vert');
title('EFG stress, \sigma_{xx}')

% Exact stresses
r=[];
theta=[];
for i=1:numnode
    sctr=node(i,:);
    Xp=sctr'; % nodal global coordinates
    xp=QT*(Xp-xTip');   % local coordinates 
    [t,rho]=cart2pol(xp(1),xp(2)); % local polar coordinates
    theta=[theta,t];
    r=[r,rho];
end
sigma=sigmato;
KI=sigma*sqrt(pi*a);


stressxx = (sigma*sqrt(pi*a)./sqrt(2*pi*r)).*cos(theta/2).*...
          (1-sin(theta/2).*sin(3*theta/2));
stressyy = (sigma*sqrt(pi*a)./sqrt(2*pi*r)).*cos(theta/2).*...
          (1+sin(theta/2).*sin(3*theta/2));
stressxy = (sigma*sqrt(pi*a)./sqrt(2*pi*r)).*sin(theta/2).*...
              cos(theta/2).*cos(3*theta/2);

figure
plot_field(node+fac*[ux' uy'],element,'Q4',stressxx);
axis('equal');
xlabel('X');
ylabel('Y');
title('\sigma_{xx}');
set(gcf,'color','white');
colorbar('vert');
title('Exact stress, \sigma_{xx}')

% ++++++++++++++++++
%    VISUALIATION
% ++++++++++++++++++

disp([num2str(toc),'   DEFORMED CONFIGURATION'])

% ----------------------------------
figure
hold on
h = plot(node(:,1)+fac*ux',...
    node(:,2)+fac*uy','r*');
set(h,'MarkerSize',7);
title('Numerical deformed shape')
axis equal
% ----------------------------------

% ----------------------------------
% Exact displacement
for i = 1 : numnode 
    x = node(i,:) ;
    [ux1,uy1] = exact_Griffith(x,E,nu,stressState,sigmato,xTip,seg,cracklength) ;
    ux_exact(i) = ux1;
    uy_exact(i) = uy1;
end
% ----------------------------------

% --------------------------------------------
% Plot both exact and numerical deformed shape
figure
hold on
h = plot(node(:,1)+fac*ux',...
    node(:,2)+fac*uy','r*');
set(h,'MarkerSize',7);
h = plot(node(:,1)+fac*ux_exact',...
    node(:,2)+fac*uy_exact','b*');
set(h,'MarkerSize',7);
title('Exact and numerical deformed shape')
legend('EFG','Exact')
axis equal
% --------------------------------------------