clear all
syms toc square_node_array  make_elem ;
D=10;
L=10;
E=1e7;
nu=0.3;
a=5;
cracklength=100;
xCr=[0 D/2;a D/2];
xTip=[a D/2];
seg=xCr(2,:)-xCr(1,:);
alpha=atan2(seg(2),seg(1));
QT=[cos(alpha) sin(alpha);-sin(alpha) cos(alpha)];
%loading
sigmato=4;

nnx=20;
nny=20;
% four corner points
pt1=[0 0];
pt2=[D 0];
pt3=[D L];
pt4=[0 L];
elemType='Q4';

numx = nnx-1;
numy = nny-1;

switch  elemType
    case 'PARTICLE' % generate the mesh for particles
       
        node=square_node_array(pt1,pt2,pt3,pt4,nnx,nny)
        inc_u=1
        inc_v=nnx
        node_pattern=[ 1 2 nnx+2 nnx+1 ]
        [element]=make_elem(node_pattern,numx,numy,inc_u,inc_v)
        
        
    case 'Q4'           % here we generate the mesh of Q4 elements
        node=square_node_array(pt1,pt2,pt3,pt4,nnx,nny)
        inc_u=1
        inc_v=nnx
        node_pattern=[ 1 2 nnx+2 nnx+1 ]
        [element]=make_elem(node_pattern,numx,numy,inc_u,inc_v)

    case 'Q9'           % here we generate a mesh of Q9 elements
        node=square_node_array(pt1,pt2,pt3,pt4,nnx,nny)
        inc_u=2
        inc_v=2*nnx
        node_pattern=[ 1 3 2*nnx+3 2*nnx+1 2 nnx+3 2*nnx+2 nnx+1 nnx+2 ]
        [element]=make_elem(node_pattern,numx,numy,inc_u,inc_v)
    case 'T3' % and last but not least T3 elements
        
        node=square_node_array(pt1,pt2,pt3,pt4,nnx,nny)
        node_pattern1=[ 1 2 nnx+1 ]
        node_pattern2=[ 2 nnx+2 nnx+1 ]
        inc_u=1
        inc_v=nnx
        element=[make_elem(node_pattern,numx,numy,inc_u,inc_v)
            make_elem(node_pattern,numx,numy,inc_u,inc_v)]
    otherwise
        error('For now, only PARTICLE, Q4, Q9 and T3 are supported by the mesh generator')
end

运行后出现Subscript indices must either be real positive integers or logicals.

Error in ==> Untitled10 at 40
        node=squarenodearray(pt1,pt2,pt3,pt4,nnx,nny)
 请高手指点 先谢过了