circular_inclusion.asv

efg code with matlab

    else
        sn = num_enr_node;
        sctrEnrB = zeros(2*(sn*1),1);
        cnt = 0 ;
        for k = 1 : le
            nodeI = index(k) ; % I-th node in index
            if (snode(k) ~= 0)
                cnt = cnt + 1 ;
                sctrEnrB(2*cnt - 1) = 2 * pos(nodeI) - 1;
                sctrEnrB(2*cnt    ) = 2 * pos(nodeI)    ;
            end
        end
        sctrB = [ sctrStdB;sctrEnrB ];
    end
   
    %  compute B matrix
    [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);

    StdB = zeros(3,2*le);
    StdB(1,1:2:2*le)  = dphidx ;
    StdB(2,2:2:2*le)  = dphidy ;
    StdB(3,1:2:2*le)  = dphidy ;
    StdB(3,2:2:2*le)  = dphidx ;

    % no enriched nodes, B is simply standard B matrix
    if (num_enr_node == 0)
        B = StdB;
    % there are enriched nodes, compute enriched part, Benr
    else
        Benr = [];
        for m = 1 : le
            if (snode(m) ~= 0)
                % compute the enrichment function and derivatives
                x = pt(1);
                y = pt(2);                
                d = sqrt((x-xc)^2+(y-yc)^2);
                F    = abs(d - r);
                dFdx = sign(F)*(x-xc)/d;
                dFdy = sign(F)*(y-yc)/d;                
                
                % B matrix of node m
                aa = dphidx(m)*F + phi(m)*dFdx ;
                bb = dphidy(m)*F + phi(m)*dFdy ;
                BI_enr = [aa 0 ; 0 bb ; bb aa];
                % Add to the total enriched B matrix
                Benr   = [Benr BI_enr];
                clear BI_enr ;
            end
        end   
        
        % Total B matrix
        B = [StdB Benr];
        clear StdB; clear Benr;
    end     
    
    % Compute the compliance matrix C
    x = pt(1);
    y = pt(2);
    d = sqrt((x-xc)^2+(y-yc)^2);
    
    if d < r 
        E  = E2 ;
        nu = nu2; 
    else
        E  = E1 ;
        nu = nu1; 
    end
    
    if ( strcmp(stressState,'PLANE_STRESS') )
        C=E/(1-nu^2)*[ 1   nu 0;
            nu  1  0 ;
            0   0  0.5*(1-nu) ];
    else
        C=E/(1+nu)/(1-2*nu)*[ 1-nu  nu  0;
            nu    1-nu 0;
            0     0  0.5-nu ];
    end
    
    % Then K matrix
    K(sctrB,sctrB) = K(sctrB,sctrB) + B'*C*B * W(igp)*J(igp) ;
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),'   INTEGRATION ON TRACTION BOUNDARY'])

% --------------------------------------
%   Generation of Gauss points
% --------------------------------------

% 4 point quadrature for each "1D element"
[W1,Q1] = quadrature(4, 'GAUSS', 1 );

Wt = [];
Qt = [];
Jt = [];

for i = 1 : (num_trac_nodes - 1)
    sctr = [tracNodes(i) tracNodes(i+1)];
    for q = 1:size(W1,1)
        pt = Q1(q,:);
        wt = W1(q);
        [N,dNdxi] = lagrange_basis('L2',pt);
        J0=dNdxi'*node(sctr,:);
        Qt = [Qt; N' * node(sctr,:)];
        Wt = [Wt; wt];
        Jt = [Jt; norm(J0)];
    end % of quadrature loop
end

% ----------------------------------
%   Form the nodal force vector f
% ----------------------------------

for igp = 1 : size(Wt,1)
    pt = Qt(igp,:);                             % quadrature point
    wt = Wt(igp);                               % quadrature weight
    detJ   = Jt(igp);
    [index] = define_support(node,pt,di);
    sctry = index.*2 ;
    [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);
    f(sctry) = f(sctry) + phi'*sigmato*detJ*wt;
end


% +++++++++++++++++++++++++++++++++++++
%   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 the bottom edge
        for i = 1 : num_disp_nodes-1      % loop on collocation points
            pt = node(dispNodes(i),:);
            % S matrix
            S = [0 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 = 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
        % for the left lower corner node, x direction
        pt = node(1,:);
        i  = i+1;
        % S matrix
        S = [1 0; 0 0]; % X 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 = 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 of if boundary point collocation is used
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,enrich_node);
    if (any(snode) == 0)
        no_enr_node  = 1 ; 
        num_enr_node = 0 ;
    else
        no_enr_node = 0;
        num_enr_node = size(find(snode ~= 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 

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

    if (no_enr_node == 1) 
        ux(i) = stdUx ;
        uy(i) = stdUy ;
    else
        % The second part exists only with enriched particles
        FUx = 0 ;
        FUy = 0 ;
        for m = 1 : le
            % a split enriched node
            if (snode(m) ~= 0)
                % compute F function
                x = pt(1);
                y = pt(2);
                d = sqrt((x-xc)^2+(y-yc)^2);
                F    = abs(d - r);
                FUx = FUx + phi(m)*F*u(2 * pos(index(m)) - 1) ;
                FUy = FUy + phi(m)*F*u(2 * pos(index(m))    ) ;
            end
        end   % end of loop on nodes in neighbour of Gp
        ux(i) = stdUx + FUx  ;
        uy(i) = stdUy + FUy ;
    end
end

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



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

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

% ----------------------------------
fac = 550 ; % 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
% ----------------------------------