pim.rar_PIM_PIM MATLAB_meshless_meshless matlab_pim method

clc % DEFINE THE PHSICAL DIMENTIONS AND MATERIAL PROPERTIES clear; length=48;height=12;ldiv=4;hdiv=4; E=30e6;nu=0.3;area=1;P=1000; quad=4;dmax=3.6; stepx=length/ldiv; stepy=height/hdiv; % DEFINE THE PLANE STRESS MATRIX C = E/(1-nu^2)*[1 nu 0;nu 1 0;0 0 (1-nu)/2]; % SET UP THE NODAL COORDINATES FOR A UNIFORM MESH [x,y,nnodes,element,elemnode]=mesh(length,height,ldiv,hdiv); [x,y]=random_func(x,y,nnodes,length,height,hdiv); %plot(x,y,'p');grid % hold on % SET UP DOMAIN OF INFLUENCE dmx=dmax*stepx*ones(1,nnodes); dmy=dmax*stepy*ones(1,nnodes); % SET UP QUADRATURE CELLS [xq,yq,numquad,qcell,elemnode]=mesh(length,height,ldiv,hdiv); % CALCULATE PQ [PQ]=PQ_func(x,y,nnodes); detPQ=det(PQ) if detPQ==0 %| detPQ==Inf | detPQ==-Inf [index,rank]=rank_func(PQ,nnodes); [PQ]=omit_func(PQ,index); detPQ=det(PQ); nnodes=rank end [cofact_PQ]=inv_func(PQ,nnodes); invPQ=cofact_PQ'/detPQ; % LOOP OVER GUASS POINTS gauss = gaussian(quad); qnodes=qcell*quad^2; [xquad,yquad,jcob,weight] = egauss(xq,yq,elemnode,gauss,qcell,nnodes); % plot(xg,yg,'p') k = zeros(2*nnodes); for j=1:qnodes xgq=xquad(j);ygq=yquad(j); [pb,dpbx,dpby]=basis_func(xgq,ygq,nnodes); [phi,dphix,dphiy]=shape_func(pb,dpbx,dpby,invPQ); BI=zeros(3,2*nnodes); for i=1:nnodes BI(1,2*i-1)=dphix(i); BI(2,2*i)=dphiy(i); BI(3,2*i-1)=dphiy(i); BI(3,2*i)=dphix(i); end k=k+BI'*C*BI*(jcob(j)*weight(j)*area); end % DETERMINE NODES ON BOUNDARY, SET UP BC'S jdis=0;jtrac=0; for i=1:nnodes if x(i)==0 jdis=jdis+1; dis_indice(jdis)=i; elseif x(i)==length jtrac=jtrac+1; trac_indice(jtrac)=i; end end %SET UP GAUSS POINTS ALONG TRACTION BOUNDARY trac_indice [xg,yg,jcob,weight] = BC_gauss(x,y,trac_indice,gauss); f=zeros(2*nnodes,1); Imo = (1/12)*height^3; qnodes=(jtrac-1)*quad; for j=1:qnodes xgq=xg(j);ygq=yg(j); [pb,dpbx,dpby]=basis_func(xgq,ygq,nnodes); [phi,dphix,dphiy]=shape_func(pb,dpbx,dpby,invPQ); if j<2 phi sum(phi) end tx=0.0; ty= -(P/(2*Imo))*((height^2)/4-ygq^2); for i=1:nnodes force(2*i-1)=phi(i)*tx; force(2*i)=phi(i)*ty; end f=f+force'*(jcob(j)*weight(j)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=1:nnodes if x(j)==0 node_temp_k=k(2*j-1,2*j-1); k(2*j-1,:)=0; k(:,2*j-1)=0; f(2*j-1,1)=0; k(2*j-1,2*j-1)=node_temp_k; end if y(j)==0 j node_temp_k=k(2*j,2*j); k(2*j,:)=0; k(:,2*j)=0; f(2*j,1)=0; k(2*j,2*j)=node_temp_k; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% u=k\f; for i=1:nnodes ux(i)=u(2*i-1); uy(i)=u(2*i); end ux uy %%%%%%%%%%%%%%%%% exact %%%%%%%%%%%%%%%%%%%% l=length*ones(size(x)); h=height*ones(size(x)); uex=P/(6*E*Imo)*((6*l-3*x).*x+(2+nu)*(y.^2-h.^2/4)) uey=-P/(6*E*Imo)*(3*nu*(y.^2).*(l-x)+(4+5*nu)*(h.^2/4).*x+(3*l-x).*(x.^2)) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55 qnodes=qcell*quad^2; [xg,yg,jcob,weight] = egauss(xq,yq,elemnode,gauss,qcell,nnodes); enorm=0; for j=1:qnodes xgq=xg(j);ygq=yg(j); [pb,dpbx,dpby]=basis_func(xgq,ygq,nnodes); [phi,dphix,dphiy]=shape_func(pb,dpbx,dpby,invPQ); BI=zeros(3,2*nnodes); for i=1:nnodes BI(1,2*i-1)=dphix(i); BI(2,2*i)=dphiy(i); BI(3,2*i-1)=dphiy(i); BI(3,2*i)=dphix(i); end stress(1:3,j) = C*BI*u; stressex(1,j) = (1/Imo)*P*(length-xgq)*ygq; stressex(2,j) = 0; stressex(3,j) = -0.5*(P/Imo)*(0.25*height^2 - ygq^2); err = stress(1:3,j)-stressex(1:3,j); err2 = weight(j)*jcob(j)*(0.5*(inv(C)*err)'*(err)); enorm = enorm + err2; end stress; enorm = sqrt(enorm)