edge_crack.m

efg code with matlab

                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 = 50 ; % visualization factor
figure
plot_field(node+fac*[ux' uy'],element,'Q4',stress_gp(1,:));
axis('equal');
title('\sigma_{xx}');
set(gcf,'color','white');
colorbar('vert');
title('EFG stress, \sigma_{xx}')
axis off

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

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

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

% -------------------------------------
%  Stress intensity factor computation
% -------------------------------------
% SIFs are computed using the interaction integral
% The integration domain is a cubic c x c, centered at
% the crack tip.
% THe weight function q is chosen q = (1-2|x|/c)(1-2|y|/c)

clear Q; clear W; clear J;

% Choose value for c; 1/4 crack length is good?
c = a/3;

% Define the integration domain
pt1 = [-c/2 -c/2] + xTip ;
pt2 = [ c/2 -c/2] + xTip ;
pt3 = [ c/2  c/2] + xTip ;
pt4 = [-c/2  c/2] + xTip ;

elemType = 'Q4' ;
[Jnode,Jelem] = meshRectangularRegion(...
    pt1, pt2, pt3, pt4,3,3,elemType);

% Build Gauss points
[w,q]=quadrature(8,'GAUSS',2);
W = [];
Q = [];
J = [];
for e = 1 : size(Jelem,1)
    sctr = Jelem(e,:);
    for i = 1:size(w,1)
        pt=q(i,:);
        wt=w(i);
        [N,dNdxi]=lagrange_basis('Q4',pt);
        J0 = Jnode(sctr,:)'*dNdxi;
        Q  = [Q; N' * Jnode(sctr,:)];
        W  = [W; wt];
        J  = [J; det(J0)];
    end
end
clear w,q ;

% -------------------------------------
% Plot the J domain
figure
hold on
cr = plot(xCr(:,1),xCr(:,2),'r-');
set(cr,'LineWidth',3)
plot_mesh(Jnode,Jelem,elemType,'b-');
plot(Q(:,1),Q(:,2),'r*');
axis equal
axis off
% -------------------------------------

% Start compute I integral

I1 = 0;
I2 = 0;
I  = [zeros(2,1)];

% -----------------------------
% start loop over Gauss points
% -----------------------------
for q = 1:size(W,1)
    pt   = Q(q,:);
    wt   = W(q);
    detJ = J(q); 
    % +++++++++++++++++++++++++
    % Gradient of displacement
    % +++++++++++++++++++++++++
    % need to compute u,x u,y v,x v,y, stored in matrix H
    [index] = define_support(node,pt,di);
    nn      = length(index);
    [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);
    leB = size(B,2);
    % nodal displacement of current element
    % taken from the total nodal parameters u
    % stdU contains true nodal displacement
    idx = 0 ;
    stdU   = zeros(2*nn,1);
    for in = 1 : nn
        idx = idx + 1;
        nodeI = index(in) ;
        stdU(2*idx-1) = u(2*nodeI-1);
        stdU(2*idx)   = u(2*nodeI  );
    end

    % A contains enriched dofs
    A = [];
    
    % number of H and tip enriched nodes
    % in index
    num_split_node = size(find(snode ~= 0),2);
    num_tip_node   = size(find(tnode ~= 0),2);

    if (num_split_node == 0) && (num_tip_node == 0)
        U = stdU ;
    else                               % having enriched DOFs
        for in = 1 : nn
            nodeI = index(in) ;
            if (snode(in) ~= 0)     % H(x) enriched node
                AA = [u(2*pos(nodeI)-1);u(2*pos(nodeI))];
                A  = [A;AA];
            elseif (tnode(in) ~= 0) % B(x) enriched node
                AA = [u(2*pos(nodeI)-1);
                    u(2*pos(nodeI));
                    u(2*(pos(nodeI)+1)-1);
                    u(2*(pos(nodeI)+1));
                    u(2*(pos(nodeI)+2)-1);
                    u(2*(pos(nodeI)+2));
                    u(2*(pos(nodeI)+3)-1);
                    u(2*(pos(nodeI)+3));
                    ];
                A  = [A;AA];
            end
        end
    end

    % total
    U = [stdU;A];

    % compute derivatives of u w.r.t xy
    H(1,1) = B(1,1:2:leB)*U(1:2:leB);    % u,x
    H(1,2) = B(2,2:2:leB)*U(1:2:leB);    % u,y
    H(2,1) = B(1,1:2:leB)*U(2:2:leB);    % v,x
    H(2,2) = B(2,2:2:leB)*U(2:2:leB);    % v,y

    % ++++++++++++++
    % Gradient of q
    % ++++++++++++++
    
    gradq(1) = -(1-2*abs(pt(2)-xTip(2))/c)*(2/c)*sign(pt(1)-xTip(1));
    gradq(2) = -(1-2*abs(pt(1)-xTip(1))/c)*(2/c)*sign(pt(2)-xTip(2));
    
    % ++++++++++++++
    % Stress at GPs
    % ++++++++++++++

    epsilon = B*U ;
    sigma   = C*epsilon;

    % +++++++++++++++++++++++++++++++++++
    % Transformation to local coordinate
    % +++++++++++++++++++++++++++++++++++

    voit2ind    = [1 3;3 2];
    gradqloc    = QT*gradq';
    graddisploc = QT*H*QT';
    stressloc   = QT*sigma(voit2ind)*QT';

    % ++++++++++++++++++
    %  Auxiliary fields
    % ++++++++++++++++++

    xp    = QT*(pt-xTip)';           % local coordinates
    r     = sqrt(xp(1)*xp(1)+xp(2)*xp(2));
    theta = atan2(xp(2),xp(1));

    K1 = 1.0 ;
    K2 = K1  ;

    mu = E/(2.+ nu + nu);
    kappa = 3-4*nu;    %Kolosov coeff, Plain strain

    SQR  = sqrt(r);
    CT   = cos(theta);
    ST   = sin(theta);
    CT2  = cos(theta/2);
    ST2  = sin(theta/2);
    C3T2 = cos(3*theta/2);
    S3T2 = sin(3*theta/2);

    drdx = CT;
    drdy = ST;
    dtdx = -ST/r;
    dtdy = CT/r;

    FACStress1 = sqrt(1/(2*pi));
    FACStress2 = FACStress1;

    FACDisp1 = sqrt(1/(2*pi))/(2*mu);
    FACDisp2 = FACDisp1;

    AuxStress   = zeros(2,2);
    AuxGradDisp = zeros(2,2);
    AuxEps      = zeros(2,2);

    for mode = 1:2
        if (mode == 1)

            AuxStress(1,1) = K1*FACStress1/SQR*CT2*(1-ST2*S3T2);
            AuxStress(2,2) = K1*FACStress1/SQR*CT2*(1+ST2*S3T2);
            AuxStress(1,2) = K1*FACStress1/SQR*ST2*CT2*C3T2;
            AuxStress(2,1) = AuxStress(1,2);

            u1    = K1*FACDisp1*SQR*CT2*(kappa - CT);
            du1dr = K1*FACDisp1*0.5/SQR*CT2*(kappa - CT);
            du1dt = K1*FACDisp1*SQR*(-0.5*ST2*(kappa - CT) + CT2*ST);

            u2    = K1*FACDisp1*SQR*ST2*(kappa - CT);
            du2dr = K1*FACDisp1*0.5/SQR*ST2*(kappa - CT);
            du2dt = K1*FACDisp1*SQR*(0.5*CT2*(kappa - CT) + ST2*ST);

            AuxGradDisp(1,1) = du1dr*drdx + du1dt*dtdx;
            AuxGradDisp(1,2) = du1dr*drdy + du1dt*dtdy;
            AuxGradDisp(2,1) = du2dr*drdx + du2dt*dtdx;
            AuxGradDisp(2,2) = du2dr*drdy + du2dt*dtdy;

            AuxEps(1,1) = AuxGradDisp(1,1);
            AuxEps(2,1) = 0.5*(AuxGradDisp(2,1) + AuxGradDisp(1,2));
            AuxEps(1,2) = AuxEps(2,1);
            AuxEps(2,2) = AuxGradDisp(2,2);
        elseif (mode == 2)
            AuxStress(1,1) = -K2*FACStress2/SQR*ST2*(2-CT2*C3T2);
            AuxStress(2,2) = K2*FACStress2/SQR*ST2*CT2*C3T2;
            AuxStress(1,2) = K2*FACStress2/SQR*CT2*(1-ST2*S3T2);
            AuxStress(2,1) = AuxStress(1,2);

            u1    = K2*FACDisp2*SQR*ST2*(kappa + 2 + CT);
            du1dr = K2*FACDisp2*0.5/SQR*ST2*(kappa + 2 + CT);
            du1dt = K2*FACDisp2*SQR*(0.5*CT2*(kappa + 2 + CT) - ST2*ST);

            u2    = -K2*FACDisp2*SQR*CT2*(kappa - 2 + CT);
            du2dr = -K2*FACDisp2*0.5*(1/SQR)*CT2*(kappa - 2 + CT);
            du2dt = -K2*FACDisp2*SQR*(-0.5*ST2*(kappa - 2 + CT) - CT2*ST);

            AuxGradDisp(1,1) = du1dr*drdx + du1dt*dtdx;
            AuxGradDisp(1,2) = du1dr*drdy + du1dt*dtdy;
            AuxGradDisp(2,1) = du2dr*drdx + du2dt*dtdx;
            AuxGradDisp(2,2) = du2dr*drdy + du2dt*dtdy;

            AuxEps(1,1) = AuxGradDisp(1,1);
            AuxEps(2,1) = 0.5*(AuxGradDisp(2,1) + AuxGradDisp(1,2));
            AuxEps(1,2) = AuxEps(2,1);
            AuxEps(2,2) = AuxGradDisp(2,2);
        end

        % +++++++++++++++
        %   J integral
        % +++++++++++++++
        I1= (stressloc(1,1) * AuxGradDisp(1,1) + stressloc(2,1) * AuxGradDisp(2,1) ) * gradqloc(1) + ...
            (stressloc(1,2) * AuxGradDisp(1,1) + stressloc(2,2) * AuxGradDisp(2,1) ) * gradqloc(2);

        I2= (AuxStress(1,1) * graddisploc(1,1) + AuxStress(2,1) * graddisploc(2,1) ) * gradqloc(1) + ...
            (AuxStress(2,1) * graddisploc(1,1) + AuxStress(2,2) * graddisploc(2,1) ) * gradqloc(2);

        StrainEnergy = 0;
        for i=1:2 
            for j=1:2  
                StrainEnergy= StrainEnergy +  stressloc(i,j)*AuxEps(i,j);
            end
        end

        % Interaction integral I
        I(mode,1) = I(mode,1) + (I1 + I2 - StrainEnergy*gradqloc(1))*detJ*wt;
    end   %loop on mode

end       % of quadrature loop
Knum = I.*E/(2*(1-nu^2)); % plain strain