edge_crack.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
% For fracture mechanics, the enriched EFG is used. This is a local PUM
% enrichment scheme which has the form:
%
%           u = N_i * u_i + N_j * H(x) * a_j + N_k * (B_i * b_ik)
%  with H(x) / Heaviside function and B_i(x) the branch functions
%  u_i, a_j and B_ik are nodal parameters.

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

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

% Problem: the infinite plate with centered crack under tension

% +++++++++++++++++++++++++++++++++++++
%          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 = 2 ;

% Material properties
E  = 1e3 ;
nu = 0.3 ;
stressState='PLANE_STRAIN';

% Loading
sigmato = 1  ;

% Crack data
a = 0.45 ;                     % crack length
xCr   = [0 D/2; a D/2];
xTip  = [a D/2];
seg   = xCr(2,:) - xCr(1,:);   % tip segment
alpha = atan2(seg(2),seg(1));  % inclination angle
QT    =[cos(alpha) sin(alpha); -sin(alpha) cos(alpha)];

% Inputs special to meshless methods, domain of influence

shape = 'circle' ;         % shape of domain of influence
dmax  = 1.7 ;              % 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
% la_approx = 100 ;       % using finite element approximation
% la_approx = 200 ;       % using point collocation method
if disp_bc_method == 1
    la_approx = 200 ;
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'])

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

% node generation, node is of size (nnode x 2)
% Corner points of the rectangle domain
% Four corner points
pt1 = [0 0] ;
pt2 = [L 0] ;
pt3 = [L D] ;
pt4 = [0 D] ;
node = square_node_array(pt1,pt2,pt3,pt4,nnx,nny);
numnode = size(node,1);

% Define boundaries
% Bottom, right and top edges are imposed by exact displacement

uln = nnx*(nny-1)+1;       % upper left node number
urn = nnx*nny;             % upper right node number
lrn = nnx;                 % lower right node number
lln = 1;                   % lower left node number
cln = nnx*(nny-1)/2+1;     % node number at (0,0)


topEdge  = [ uln:1:(urn-1); (uln+1):1:urn ]';
botEdge  = [ lln:1:(lrn-1); (lln+1):1:lrn ]';

% GET NODES ON DIRICHLET BOUNDARY AND ESSENTIAL BOUNDARY
botNodes  = unique(botEdge);
topNodes  = unique(topEdge);

dispNodes = [botNodes];
tracNodes = topNodes;

% The bottom edge is constrained along Y direction
% The first node of this edge is constrained to X direction
num_disp_nodes = length(dispNodes)+1;
num_trac_nodes = length(tracNodes)  ;

% 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

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
[w,q]=quadrature(6, 'GAUSS', 2 );  % 4x4 Gaussian quadrature for each cell
W = [];
Q = [];
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 ;

% +++++++++++++++++++++++++++++++++++++
%        LEVEL SET INITIALIZATION
%      SELECTION OF ENRICHED NODES
% +++++++++++++++++++++++++++++++++++++
% Only normal level sets are computed and stored in nodes
x0  = xCr(1,1); y0 = xCr(1,2);
x1  = xCr(2,1); y1 = xCr(2,2);
t   = 1/norm(seg)*seg;
for i = 1 : numnode
    x = node(i,1);
    y = node(i,2);
    l   = sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)) ;
    phi = (y0-y1)*x + (x1-x0)*y + (x0*y1-x1*y0);
    ls(i,1) = phi/l;            % normal LS
    ls(i,2) = ([x y]-xTip)*t';  % tangent LS
    dt(i)   = norm([x y]-xTip); % distances from node I to tip
end

set1 = find(abs(ls(:,1)) - di' < 0);
set2 = find(ls(:,2) < 0);
set3 = find(dt - di < 0);        % tip nodes
set4 = intersect(set1,set2);     % split nodes

% some nodes belong to both sets, remove tip
split_nodes = setdiff(set4,set3);
tip_nodes   = set3;

% +++++++++++++++++++++++++++++++++++++
%  PLOT NODES,BACKGROUND MESH, GPOINTS
% +++++++++++++++++++++++++++++++++++++
disp([num2str(toc),'   PLOT NODES AND GAUSS POINTS'])
figure
hold on
cntr = plot([0,L,L,0,0],[0,0,D,D,0]);
set(cntr,'LineWidth',3)
h = plot(node(:,1),node(:,2),'bo');
set(h,'MarkerSize',8.5);
cr = plot(xCr(:,1),xCr(:,2),'r-');
set(cr,'LineWidth',3)
n1 = plot(node(split_nodes,1),node(split_nodes,2),'r*');
n2 = plot(node(tip_nodes,1),node(tip_nodes,2),'rh');
set(n1,'MarkerSize',15);
set(n2,'MarkerSize',15);
plot(node(dispNodes,1),node(dispNodes,2),'ks');
axis off

%theta = -pi:0.1:pi;
%for i = 1 : size(split_nodes',1)
%    a = di(i); % radius of support
%    x = node(split_nodes(i),1) + a*cos(theta);
%    y = node(split_nodes(i),2) + a*sin(theta);
%    cir = plot(x,y,'k-');
%end

%for i = 1 : size(set3',1)
%    a = di(i); % radius of support
%    x = node(set3(i),1) + a*cos(theta);
%    y = node(set3(i),2) + a*sin(theta);
%    cir = plot(x,y,'r.-');
%end

figure
hold on
cntr = plot([0,L,L,0,0],[0,0,D,D,0]);
set(cntr,'LineWidth',3)
plot_mesh(node,element,'Q4','--');
plot(Q(:,1),Q(:,2),'r*');
cr = plot(xCr(:,1),xCr(:,2),'r-');
set(cr,'LineWidth',3)
axis off

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

% Initialisation

total_num_node = numnode + size(split_nodes,2)*1 + size(tip_nodes,2)*4 ;
K = sparse(2*total_num_node,2*total_num_node);
f = zeros(2*total_num_node,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

% Due to the presence of additional dofs, the assembly is a little
% bit difficult than in FEM. We use fictitious nodes to handle these
% additional dofs. At a H(x) enriched node, we add one fantom node and
% at tip enriched node, four fantom nodes are added. These fictitious nodes
% are numbered from the total number of true nodes, ie, from numnode+1 ...

pos = zeros(numnode,1);
nsnode = 0 ;
ntnode = 0 ;
for i = 1 : numnode
    [snode] = ismember(i,split_nodes);
    [tnode] = ismember(i,tip_nodes);
    if (snode ~= 0)
        pos(i) = (numnode + nsnode*1 + ntnode*4) + 1 ;
        nsnode = nsnode + 1 ;
    elseif (tnode ~= 0)
        pos(i) = (numnode + nsnode*1 + ntnode*4) + 1 ;
        ntnode = ntnode + 1 ;
    end
end


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

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

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

    % For enriched EFG, here is the most important part
    % Nodes within neighbouring of GPs include normal nodes, H enriched
    % nodes and tip enriched nodes
    % ISMEMBER(A,S) for the array A returns an array of the same size as A
    % containing 1 where the elements of A are in the set S and 0 otherwise.

    [sctrB,snode,tnode] = assembly(index,split_nodes,tip_nodes,...
        pos);
    %  compute B matrix
    B = Bmatrix(pt,index,node,di,form,...
        snode,tnode,xCr,xTip,alpha);

    % 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


    % if Lagrange multipliers are interpolated by finite elements
    % along the displacement boundary
    if (la_approx == 100)
        disp(' Finite element interpolation for ML is used ')

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

        % ------------------------------
        %   Generation of Gauss points
        % ------------------------------
        Wu = [];
        Qu = [];
        Ju = [];

        for i = 1 : (num_disp_nodes - 1)
            sctr = [disp_nodes(i) disp_nodes(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];