plate_hole.m

efg code with matlab

clear Wt; clear Qt; clear Jt;

% +++++++++++++++++++++++++++++++++++++
%   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

    % ---------------------------------
    %    Form vector q and matrix G
    % ---------------------------------
    qk = zeros(1,2*num_disp_nodes);
    G  = zeros(2*numnode,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 ')
        ind = 0 ;
        for i = 1 : ndnode_1      % loop on collocation points
            ind = ind + 1;
            pt = node(d_node_1(i),:);
            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*ind-1:2*ind) = G(row1:row2,2*ind-1:2*ind) + Gj ;
            end
        end

        % The same thing for the left edge

        for i = 1 : ndnode_2      % loop on collocation points
            ind = ind + 1;
            pt = node(d_node_2(i),:);
            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*ind-1:2*ind) = G(row1:row2,2*ind-1:2*ind) + 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 : (ndnode_1 - 1)
            sctr = [d_node_1(i) d_node_1(i+1)];
            m1 = m1 + 1 ;
            m2 = m1 + 1 ;
            le = norm(d_node_1(m1) - d_node_1(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
                N1 = 1 - pt(1,1)/le ; N2 = 1 - N1 ;
                % 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)   ;
                    G1   = - wt * detJ * phi(j) * N1 * S ;
                    G2   = - wt * detJ * phi(j) * N2 * S ;
                    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 of loop on Gauss points
        end  % end of loop on 1D elements

        % The same thing for the left edge
        m1 = 0 ;
        for i = 1 : (ndnode_2 - 1)
            sctr = [d_node_2(i) d_node_2(i+1)];
            m1 = m1 + 1 ;
            m2 = m1 + 1 ;
            le = norm(d_node_2(m1) - d_node_2(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
                N1 = 1 - pt(1,1)/le ; N2 = 1 - N1 ;
                % 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)   ;
                    G1   = - wt * detJ * phi(j) * N1 * S ;
                    G2   = - wt * detJ * phi(j) * N2 * S ;
                    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 of loop on Gauss points
        end  % end of loop on 1D elements

    end   % end of finite element interpolation for ML

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 % end of if Penalty method

% +++++++++++++++++++++++++++++++++++++
%      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}
    u2 = d(1:numnode*2);
    % just get displacement parameters u_i
    for i = 1 : numnode
        u(1,i) = d(2*i-1); % x displacement
        u(2,i) = d(2*i);   % y displacement
    end
end

clear d ;

% +++++++++++++++++++++++++++++++++++++
%    COMPUTE THE TRUE DISPLACEMENTS
% +++++++++++++++++++++++++++++++++++++
disp([num2str(toc),'   INTERPOLATION TO GET TRUE DISPLACEMENT'])
disp = zeros(1,2*numnode);
for i = 1 : numnode
    [index] = define_support(node,node(i,:),di);
    % shape function at nodes in neighbouring of node i
    [phi,dphidx,dphidy] = MLS_ShapeFunction(node(i,:),index,node,di,form);
    disp(1,2*i-1) = phi*u(1,index)'; % x nodal displacement
    disp(1,2*i)   = phi*u(2,index)'; % y nodal displacement
end
clear u;
% +++++++++++++++++++++++++++++++++++++
%   COMPUTE STRESSES AT GAUSS POINTS
% +++++++++++++++++++++++++++++++++++++
disp([num2str(toc),'   COMPUTE STRESS AT GAUSS POINTS'])
ind = 0 ;
for igp = 1 : size(W,1)
    pt = Q(igp,:);
    [index] = define_support(node,pt,di);
    B = zeros(3,2*size(index,2)) ;
    en = zeros(1,2*size(index,2));
    [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);
    for m = 1 : size(index,2)
        B(1:3,2*m-1:2*m) = [dphidx(m) 0 ;
            0 dphidy(m);
            dphidy(m) dphidx(m)];
        en(2*m-1) = 2*index(m)-1;
        en(2*m  ) = 2*index(m)  ;
    end
    ind = ind + 1 ;
    stress_gp(1:3,ind) = C*B*u2(en);  % sigma = C*epsilon = C*(B*u)
end

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

% Deformed configuration
% ----------------------

figure
plot_mesh(node,element,'Q4','k--');
fac = 100 ; % visualization factor
hold on
h = plot(node(:,1)+fac*disp(1,1:2:2*numnode-1)',...
    node(:,2)+fac*disp(1,2:2:2*numnode)','r*');
set(h,'MarkerSize',7);
axis off
title('Deformed configuration (scaled)');

% Stress visualization
% ----------------------
% Stresses are computed at Gauss points
% Gauss points are then used to build a Delaunay triangulation, then plot
figure
hold on
tri = delaunay(Q(:,1),Q(:,2));
plot_field(Q,tri,'T3',stress_gp(1,:));
axis('equal');
xlabel('X');
ylabel('Y');
title('Sigma XX');
set(gcf,'color','white');
colorbar('vert');
opts = struct('Color','rgb','Bounds','tight');
exportfig(gcf,'beam_stress_x.eps',opts)

% Export to Tecplot for better plot
varname{1}='X';
varname{2}='Y';
varname{3}='sigma_xx';
varname{4}='sigma_yy';
varname{5}='sigma_xy';
%varname{6}='sigma_vonMises';
sigma=[stress_gp(1,:)' stress_gp(2,:)' stress_gp(3,:)'];
value=[Q(:,1) Q(:,2) sigma];
tecplotout(tri,'T3',Q,varname,value,'hole_stress.dat');

for igp = 1 : size(W,1)
    pt = Q(igp,:);
    str = exact_plate_hole(pt,a);
    exact_sigma(:,igp) = str';
end

varname{1}='X';
varname{2}='Y';
varname{3}='sigma_xx';
value=[Q(:,1) Q(:,2) sigma_xx'];
tecplotout(tri,'T3',Q,varname,value,'beam_exact_stress.dat');

% Remove used memory
clear coord; clear Q; clear W; clear J; clear stress_gp;

% --------------------------------------------------------------------
%                           END OF THE PROGRAM
% --------------------------------------------------------------------
figure
plot_field(Q,tri,'T3',exact_sigma(1,:));
axis('equal');
xlabel('X');
ylabel('Y');
title('Sigma XX');
set(gcf,'color','white');
colorbar('vert');

figure
plot_field(Q,tri,'T3',exact_sigma(2,:));
axis('equal');
xlabel('X');
ylabel('Y');
title('Sigma XX');
set(gcf,'color','white');
colorbar('vert');


ind = 0 ;
for igp = 1 : size(node,1)
    pt = node(igp,:);
    wt = W(igp);
    [index] = define_support(node,pt,di);
    B  = zeros(3,2*size(index,2)) ;
    en = zeros(1,2*size(index,2));
    [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);
    for m = 1 : size(index,2)
        B(1:3,2*m-1:2*m) = [dphidx(m) 0 ;
            0 dphidy(m);
            dphidy(m) dphidx(m)];
        en(2*m-1) = 2*index(m)-1;
        en(2*m  ) = 2*index(m)  ;
    end
    ind = ind + 1 ;
    stress_gp(1:3,ind) = C*B*u2(en);  % sigma = C*epsilon = C*(B*u)
end

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

for igp = 1 : size(node,1)
    pt = node(igp,:);
    str = exact_plate_hole(pt,a);
    exact_sigma(:,igp) = str';
end

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