problem_boundary_condition.m

efg code with matlab

        m1 = m1 + 1 ;
        m2 = m1 + 1 ;
        le = norm(node(sctr(1),:) - node(sctr(2),:));
        %le = norm( disp_nodes_3(i,:) - disp_nodes_3(i+1,:) );
        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
            %pt =  N' * disp_nodes_3(sctr,:); % global GP
            N1 = 1 - pt(1,1)/le ; N2 = 1 - N1 ;
            % qk vector
            qk(2*m1-1) = qk(2*m1-1) - wt * detJ * N1 * 0 ;
            qk(2*m1)   = qk(2*m1)   - wt * detJ * N1 * (-0.2) ;
            qk(2*m2-1) = qk(2*m2-1) - wt * detJ * N2 * 0 ;
            qk(2*m2)   = qk(2*m2)   - wt * detJ * N2 * (-0.2) ;
            % G matrix
            [index] = define_support(node,pt,di);
            [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);
            % S matrix, BCs imposed on y direction
            S = [0 0 ; 0 1];
            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 GPs of each element
    end % end of loop on 1D elements
end
end

% *****************************
%   Penalty function is used
% *****************************

% 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

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

if disp_bc_method == 2
    fu = zeros(2*numnode,1);
    k = zeros(2*numnode,2*numnode);
    % for bottom edge
    
    % ------------------------------
    %   Generation of Gauss points
    % ------------------------------
    Wu = [];
    Qu = [];
    Ju = [];

    for i = 1 : (size(disp_nodes_1,1) - 1)
        sctr = [disp_nodes_1(i) disp_nodes_1(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
    
    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);
        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)*0;
            force(2*j  ) = phi(j)*0;
        end
        fu(en) = fu(en) + Ju(igp) * wt * force'   ;
        phi_ij = phi_ij + Ju(igp) * wt * phi*phi' ;
    end % end of loop on Gauss points
    f = f - alpha*fu;
    K(en,en) = K(en,en) - alpha*phi_ij ; 
    
    % for left edge
    
    % ------------------------------
    %   Generation of Gauss points
    % ------------------------------
    Wu = [];
    Qu = [];
    Ju = [];
    
    for i = 1 : (size(disp_nodes_2,1) - 1)
        sctr = [disp_nodes_2(i) disp_nodes_2(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
    
    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);
        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)*0;
            force(2*j  ) = phi(j)*0;
        end
        fu(en) = fu(en) + Ju(igp) * wt * force' ;
        phi_ij = phi_ij + Ju(igp)*wt*phi*phi' ;
    end % end of loop on Gauss points
    f = f - alpha*fu;
    K(en,en) = K(en,en) - alpha*phi_ij ; 
    
    % for top edges
    
    % ------------------------------
    %   Generation of Gauss points
    % ------------------------------
    Wu = [];
    Qu = [];
    Ju = [];

    for i = 1 : (size(disp_nodes_3,1) - 1)
        sctr = [disp_nodes_3(i) disp_nodes_3(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
   
    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);
        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)*0;
            force(2*j  ) = phi(j)*(-0.2);
        end
        fu(en) = fu(en) + Ju(igp) * wt * force' ;
        phi_ij = phi_ij + Ju(igp)*wt*phi*phi' ;
    end % end of loop on Gauss points
    
    % do the subtract
    f = f - alpha*fu;
    K(en,en) = K(en,en) - alpha*phi_ij ; 
    
end

% +++++++++++++++++++++++++
%   Modified variational
% +++++++++++++++++++++++++

if disp_bc_method == 3
    fu = zeros(2*numnode,1);
    k = zeros(2*numnode,1);
    % for bottom edge
    
    % ------------------------------
    %   Generation of Gauss points
    % ------------------------------
    Wu = [];
    Qu = [];
    Ju = [];

    for i = 1 : (size(disp_nodes_1,1) - 1)
        sctr = [disp_nodes_1(i) disp_nodes_1(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
    
    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);
        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)*0;
            force(2*j  ) = phi(j)*0;
        end
        fu(en) = fu(en) + Ju(igp) * wt * force'   ;
        phi_ij = phi_ij + Ju(igp) * wt * phi*phi' ;
    end % end of loop on Gauss points
    f = f - alpha*fu;
    K(en,en) = K(en,en) - alpha*phi_ij ; 
    
    % for left edge
    
    % ------------------------------
    %   Generation of Gauss points
    % ------------------------------
    Wu = [];
    Qu = [];
    Ju = [];
    
    for i = 1 : (size(disp_nodes_2,1) - 1)
        sctr = [disp_nodes_2(i) disp_nodes_2(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
    
    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);
        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)*0;
            force(2*j  ) = phi(j)*0;
        end
        fu(en) = fu(en) + Ju(igp) * wt * force' ;
        phi_ij = phi_ij + Ju(igp)*wt*phi*phi' ;
    end % end of loop on Gauss points
    f = f - alpha*fu;
    K(en,en) = K(en,en) - alpha*phi_ij ; 
    
    % for top edges
    
    % ------------------------------
    %   Generation of Gauss points
    % ------------------------------
    Wu = [];
    Qu = [];
    Ju = [];

    for i = 1 : (size(disp_nodes_3,1) - 1)
        sctr = [disp_nodes_3(i) disp_nodes_3(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
   
    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);
        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)*0;
            force(2*j  ) = phi(j)*(-0.2);
        end
        fu(en) = fu(en) + Ju(igp) * wt * force' ;
        phi_ij = phi_ij + Ju(igp)*wt*phi*phi' ;
    end % end of loop on Gauss points
    
    % do the subtract
    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 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

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

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;

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

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

figure
plot_mesh(node,element,'Q4','k--');
fac = 1 ; % 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);
title('Deformed configuration (scaled)');
opts = struct('Color','rgb','Bounds','tight');
exportfig(gcf,'boundary_deformed3.eps',opts);


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

% ------------------------------
%      END OF THE PROGRAM
% ------------------------------