plate_hole.asv

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

% 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  = 4;                 % 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 = 1 ;

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

% If penalty function is used, then choose "smart" penalty number :-)
if disp_bc_method == 2
    alpha = 1e7 ;           % 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'])

% Name of the mesh file
meshFile = 'hole.asc';
% Read this file to get nodes and even elements as integration cells
[node,element,elemType] = sysmesh2mlab(meshFile);
numnode = size(node,1);
node = node(:,1:2);  % two dimension, no z coordinate

% Boundary nodes
% 1. Displacement conditions:
%    - On bottom edge (y=0) with u_y = 0
%    - On left   edge (x=0) with u_x = 0

% 2. Traction conditions
%    - On right edge (x=5) with t = (sigma_x,sigma_xy)
%    - On top edge   (y=5) with t = (sigma_xy,sigma_y)

% from the file hole.asc
% Need to write the code to read instead of copy like this :-)
d_node_1 = [11 10 9 8 7 6 5 4 3 2 1] ;
d_node_2 = [99 98 97 96 95 94 93 92 91 90 89] ;
ndnode_1 = length(d_node_1);
ndnode_2 = length(d_node_2);
num_disp_nodes = ndnode_1 + ndnode_2;

t_node_1 = [45 34 23 12 1] ;
t_node_2 = [89 78 67 56 45] ;
ntnode_1 = length(t_node_1);
ntnode_2 = length(t_node_2);

%ubar = zeros(2*num_disp_nodes,1); % initialize of vector u bar
f    = zeros(2*numnode,1);        % nodal force vector

% Domain of influence for every nodes 
di = zeros(1,numnode);
for e = 1 : size(element,1)        % element loop
    sctr=element(e,:);             % element scatter vector
    % compute the diagonal of this subcell e
    node1 = node(sctr(1),:);
    node3 = node(sctr(3),:);
    d = norm(node1-node3);
    di(1,sctr) = d*dmax;
end

% +++++++++++++++++++++++++++++++++++++
%            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 '])

% 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(elemType,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'])
% plot nodes
figure
hold on
h = plot(node(:,1),node(:,2),'bo');
set(h,'MarkerSize',9);
cd = plot([1 5 5 0 0],[0 0 5 5 1],'b-');
set(cd,'LineWidth',2);
theta = 0:0.1:pi/2;
a = 1; % radius of hole
x = a*cos(theta);
y = a*sin(theta);
cir = plot(x,y,'b-');
set(cir,'LineWidth',2);
plot(node(d_node_1,1),node(d_node_1,2),'r*');
plot(node(d_node_2,1),node(d_node_2,2),'r*');
plot(node(t_node_1,1),node(t_node_1,2),'bs');
plot(node(t_node_2,1),node(t_node_2,2),'bs');
axis equal
axis off
axis([-1 6 -1 6])
title('Nodes distribution');
opts = struct('Color','rgb','Bounds','tight');
exportfig(gcf,'hole_nodes.eps',opts);

% plot background mesh and Gauss points
figure
hold on
cd = plot([1 5 5 0 0],[0 0 5 5 1],'b-');
set(cd,'LineWidth',2);
plot_mesh(node,element,'Q4','--');
plot(Q(:,1),Q(:,2),'r*');
axis([-1 6 -1 6])
title('Gauss points');
axis off
exportfig(gcf,'hole_gp.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 TRACTION BOUNDARYS
% +++++++++++++++++++++++++++++++++++++
% 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'])

% Since we have two traction boundaries, we do two times the procedure

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

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

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

for i = 1 : (ntnode_1 - 1)
    sctr = [t_node_1(i) t_node_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,:);                    
        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
    [index] = define_support(node,pt,di);
    en = zeros(1,2*size(index,2));              
    force = zeros(1,2*size(index,2));
    [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);
    
    % compute the exact solution 
    str = exact_plate_hole(pt,a);
    % prescribed traction vector t_bar
    tx = str(1) ;
    ty = str(3);
    for j = 1 : size(index,2)
        en(2*j-1) = 2*index(j)-1;
        en(2*j  ) = 2*index(j)  ;
        force(2*j-1) = tx*phi(j);
        force(2*j  ) = ty*phi(j);
    end
    f(en) = f(en) + Jt(igp) * wt * force' ;
end

clear Wt; clear Qt; clear Jt; 

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

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

for i = 1 : (ntnode_2 - 1)
    sctr = [t_node_2(i) t_node_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,:);                    
        Qt = [Qt; N' * node(sctr,:)];              
        Wt = [Wt; wt];                             
        Jt = [Jt; norm(J0)];
    end % of quadrature loop
end