problem_boundary_condition.m

efg code with matlab

% *************************************************************************
%                 TWO DIMENSIONAL ELEMENT FREE GALERKIN CODE
%                            Nguyen Vinh Phu
%                        LTDS, ENISE, Juillet 2006
% *************************************************************************

% Description:
% This is a simple Matlab code of the EFG method which is the most familiar
% meshless methods using MLS approximation and Galerkin procedure to derive
% the discrete equations.
% Domain of influence can be : (1) circle and (2) rectangle
% Weight function : cubic or quartic spline function
% Nodes can be uniformly distributed nodes or randomly distributed
% (in the latter case, this is read from a finite element mesh file)
% Numerical integration is done with background mesh
% Essential boundary condition is imposed using Lagrange multipliers or
% penalty method.
% In case of Lagrange multipliers, but use the Dirac delta function to
% approximate lamda, we obtain the boundary point collocation method

% External utilities: the mesh, post processing is done with the Matlab
% code of Northwestern university, Illinois, USA

% *************************************************************************

% Problem: Infinite plate with centered hole
% Plate dimensions: 10x10 with hole of unity radius
% The plate is meshed with SYSTUS program and stored in the file hole.asc
% Due to symmetry, only 1/4 plate is modeled.

% +++++++++++++++++++++++++++++++++++++
%          PROBLEM INPUT DATA
% +++++++++++++++++++++++++++++++++++++

clear all
clc
state = 0;
tic;           % help us to see the time required for each step

% Dimension of the domain (it is simply a rectangular region L x W)
L = 1 ;
D = 1 ;

% Material properties
E  = 1000;
nu = 0.3 ;
stressState='PLANE_STRAIN'; % set to either 'PLANE_STRAIN' or "PLANE_STRESS'

% Inputs special to meshless methods, domain of influence

shape = 'circle' ;         % shape of domain of influence
dmax  = 1.5;                 % radius = dmax * nodal spacing
form  = 'cubic_spline' ;   % using cubic spline weight function

% Choose method to impose essential boundary condition
% disp_bc_method = 1 : Lagrange multiplier method
% disp_bc_method = 2 : Penalty method
% disp_bc_method = 3 : Modified variational principles

disp_bc_method = 1 ;

% If Lagrange multiplier is used, choose approximation to be used for lamda
% la_approx = 100 ;       % using finite element approximation
% la_approx = 200 ;       % using point collocation method   
if disp_bc_method == 1
    la_approx = 100 ;         
end

% If penalty function is used, then choose "smart" penalty number :-)
if disp_bc_method == 2
    alpha = 1e5 ;           % penalty number
end

% +++++++++++++++++++++++++++++++++++++
%            NODE GENERATION
% +++++++++++++++++++++++++++++++++++++
% Steps:
%  1. Node coordinates
%  2. Nodes' weight function including : shape (circle or rectangle), size
%  and form (quartic spline or the other)

disp([num2str(toc),'   NODE GENERATION'])

% Node density defined by number of nodes along two directions
nnx = 13 ;
nny = 13  ;

% node generation, node is of size (nnode x 2)
node = square_node_array([0 0],[1 0],[1 1],[0 1],nnx,nny);
numnode = size(node,1);

% Boundary nodes
% 1. Displacement condition on the left edge with x = 0
% 2. Traction condition on the right edge with x = L

disp_nodes_1 = find(node(:,2)==0) ;
disp_nodes_2 = find(node(:,1)==0) ;
disp_nodes_3 = intersect(find(node(:,2)==L),find(node(:,1) <= L/3) ) ;


num_disp_nodes = length(disp_nodes_1)+length(disp_nodes_2)+...
                 length(disp_nodes_3);
     
% define collocation points along displacement boundaries

%disp_nodes_1 = linspace(0,D,24)';
%disp_nodes_1 = [disp_nodes_1 zeros(length(disp_nodes_1),1)];

%disp_nodes_2 = linspace(0,L,24)';
%disp_nodes_2 = [zeros(length(disp_nodes_2),1) disp_nodes_2 ];

%disp_nodes_3 = linspace(0,L/3,8)';
%disp_nodes_3 = [disp_nodes_3 ones(length(disp_nodes_3),1)];

%num_disp_nodes = size(disp_nodes_1,1)+size(disp_nodes_2,1)+...
%                 size(disp_nodes_3,1);

% Domain of influence for every nodes 
% Uniformly distributed nodes
% Definition : rad = dmax*max(deltaX,deltaY)
% deltaX,deltaY are nodal spacing along x,y directions

deltaX = L/(nnx-1);
deltaY = D/(nny-1);
delta  = max(deltaX,deltaY);
di     = ones(1,numnode)*dmax*delta ;

% +++++++++++++++++++++++++++++++++++++
%            GAUSS POINTS
% +++++++++++++++++++++++++++++++++++++
% Steps:
% 1. Mesh generation using available meshing algorithm of FEM
% 2. Using isoparametric mapping to transform Gauss points defined in each
%   element (background element) to global coordinates
% For this problem, since the geometry is not regular, the background mesh
% is used.

disp([num2str(toc),'   BUILD GAUSS POINT FOR DOMAIN '])

% Meshing, Q4 elements
inc_u = 1;
inc_v = nnx;
node_pattern = [ 1 2 nnx+2 nnx+1 ];
element = make_elem(node_pattern,nnx-1,nny-1,inc_u,inc_v);

% Building Gauss points
% 4x4 Gaussian quadrature for each element
[w,q]=quadrature(4, 'GAUSS', 2 );
W = [];   % weight
Q = [];   % coord
J = [];   % det J
for e = 1 : size(element,1)        % element loop
    sctr=element(e,:);             % element scatter vector
    for i = 1:size(w,1)                        % quadrature loop
        pt=q(i,:);                             % quadrature point
        wt=w(i);                               % quadrature weight
        [N,dNdxi]=lagrange_basis('Q4',pt);     % Q4 shape functions
        J0=node(sctr,:)'*dNdxi;                % element Jacobian matrix
        Q = [Q; N' * node(sctr,:)];            % global Gauss point
        W = [W; wt];                           % global Gauss point
        J = [J; det(J0)];
    end
end
clear w,q ;

% +++++++++++++++++++++++++++++++++++++
%  PLOT NODES,BACKGROUND MESH, GPOINTS
% +++++++++++++++++++++++++++++++++++++
disp([num2str(toc),'   PLOT NODES AND GAUSS POINTS'])

disp_nodes = [disp_nodes_1;disp_nodes_2;disp_nodes_3];
%disp_nodes = [disp_nodes_1;disp_nodes_2];

figure
hold on
h = plot(node(:,1),node(:,2),'bo');
set(h,'MarkerSize',6);
cd = plot([0 1 1 0 0],[0 0 1 1 0],'b-');
set(cd,'LineWidth',2);
g = plot(node(disp_nodes,1),node(disp_nodes,2),'rp');
set(g,'MarkerSize',12);
%s = plot(disp_nodes_1(:,1),disp_nodes_1(:,2),'rp');
%set(s,'MarkerSize',12);
%s = plot(disp_nodes_2(:,1),disp_nodes_2(:,2),'rp');
%set(s,'MarkerSize',12);
%s = plot(disp_nodes_3(:,1),disp_nodes_3(:,2),'rp');
%set(s,'MarkerSize',12);
axis equal
axis off

opts = struct('Color','rgb','Bounds','tight');
exportfig(gcf,'boundary_collocation_point2.eps',opts);

% +++++++++++++++++++++++++++++++++++++
%          DOMAIN ASSEMBLY
% +++++++++++++++++++++++++++++++++++++

disp([num2str(toc),'   DOMAIN ASSEMBLY'])

% Initialisation
K = sparse(2*numnode,2*numnode);
f = zeros(2*numnode,1);

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


% ---------------------------------
%       Loop on Gauss points
% ---------------------------------

for igp = 1 : size(W,1)
    pt = Q(igp,:);               % quadrature point
    wt = W(igp);                 % quadrature weight

    % ----------------------------------------------
    % find nodes in neighbouring of Gauss point pt
    % ----------------------------------------------
    [index] = define_support(node,pt,di);

    % --------------------------------------
    %   compute B matrix, K matrix
    % --------------------------------------
    % Loop on nodes within the support, compute dPhi of each node
    % Then get the B matrix
    % Finally assemble to the K matrix

    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
    K(en,en) = K(en,en) + B'*C*B * W(igp)*J(igp) ;
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 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);

% *****************************
%   Point collocation is used
% *****************************

if (la_approx == 200)

    % bottom edge with imposed BCs on u_y
    % +++++++++++++++++++++++++++++++++++
    ind = 0 ;
    for i = 1 : size(disp_nodes_1,1)
        ind = ind + 1 ;
        pt  = disp_nodes_1(i,:);
        S   = [0 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*ind-1:2*ind) = G(row1:row2,2*ind-1:2*ind) + Gj ;
        end
    end

    % left edge with imposed BCs on u_x
    % +++++++++++++++++++++++++++++++++
    for i = 1 : size(disp_nodes_2,1)
        ind = ind + 1 ;
        pt  = disp_nodes_2(i,:);
        S   = [1 0; 0 0]; % 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*ind-1:2*ind) = G(row1:row2,2*ind-1:2*ind) + Gj ;
        end
    end

    % top edge with imposed BCs on u_y = - 0.2
    % ++++++++++++++++++++++++++++++++++++++++
    for i = 1 : size(disp_nodes_3,1)
        ind = ind + 1 ;
        pt = disp_nodes_3(i,:);
        S = [0 0; 0 1]; % both directions
        % qk vector
        qk(1,2*ind-1) = 0 ;
        qk(1,2*ind)   = 0.2 ;
        % 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

% *****************************
%   FE interpolation for ML
% *****************************

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

if (la_approx == 100)
    disp(' Finite element interpolation for ML is used ')
    m1 = 0 ;
    for i = 1 : (size(disp_nodes_1,1) - 1)
        sctr = [disp_nodes_1(i) disp_nodes_1(i+1)];
        m1 = m1 + 1 ;
        m2 = m1 + 1 ;
        le = norm(node(sctr(1),:) - node(sctr(2),:));
        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 ;
            % 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
    
    % For the left edge ...
    for i = 1 : (size(disp_nodes_2,1) - 1)
        sctr = [disp_nodes_2(i) disp_nodes_2(i+1)];
        m1 = m1 + 1 ;
        m2 = m1 + 1 ;
        le = norm(node(sctr(1),:) - node(sctr(2),:));
        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,2)/le ; N2 = 1 - N1 ;
            % G matrix
            [index] = define_support(node,pt,di);
            [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);
            % S matrix, BCs imposed on x direction
            S = [1 0 ; 0 0];
            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
    
    % for the top edge
    
    for i = 1 : (size(disp_nodes_3,1) - 1)
        %sctr = [i i+1];
        sctr = [disp_nodes_3(i) disp_nodes_3(i+1)];