MLPG Patch.rar_MESHFREE_MESHFREE METHOD_MLPG

% PROGRAM FOR THE MESHLESS RIGID BODY FORMULATION % 2D - PATCH TEST clc; clear; % DEFINE BOUNDARIES/PARAMETERS Lb = 1; D = 1; young = 1;nu=0.25; P=1; tor=1e-10; ng=4; %number of gauss points m=3; % PLANE STRESS MATRIX Dmat = (young/(1-nu^2))*[1 nu 0;nu 1 0;0 0 (1-nu)/2]; % SET UP NODAL COORDINATES ndivl=2; ndivw=2; % [Coords,conn,numcell,numnod,id] = mesh2d(Lb,D,ndivl,ndivw); Coords=[0 0 0 .5 .5 .5 1 1 1;1 .5 0 1 .5 0 1 .5 0] % Coords=[0 0 0 .666667 .5 .333333 1 1 1;1 .6666667 0 1 .5 0 1 .333333 0] dlmwrite('coords_r2.dat',Coords); numnod=size(Coords,2); % DETERMINE DOMAINS OF INFLUENCE - UNIFORM NODAL SPACING alfs=4.5; alfq=0.5; xspac = Lb/ndivl; yspac = D/ndivw; ds(1,1:numnod)=alfs*xspac*ones(1,numnod); ds(2,1:numnod)=alfs*yspac*ones(1,numnod); % ESSENTIAL BOUNDARY CONDITIONS ebc(1,1)=3; %Node 3 ebc(1,2)=1; % Flag x direction ebc(1,3)=1; % Flag y direction ebc(1,4)=0; ebc(1,5)=0; ebc(2,1)=6; %Node 6 ebc(2,2)=0; ebc(2,3)=1; ebc(2,4)=0; ebc(2,5)=0; ebc(3,1)=9; %Node 9 ebc(3,2)=0; ebc(3,3)=1; ebc(3,4)=0; ebc(3,5)=0; %LOOP ALL FIELD NODES - ASSEMBLY SYSTEM OF EQUATIONS kg=zeros(2*numnod,2*numnod); fg=zeros(1,2*numnod); for nod=1:numnod lin1=(nod*2)-1; lin2=nod*2; QDCoords = Qdomain( alfq*xspac, alfq*yspac, Coords(1,nod), Coords(2,nod), Coords); for seg=1:4 % Loop over quad. domain segments if (seg+1)==5 xi=[QDCoords(1,seg);QDCoords(2,seg)]; xj=[QDCoords(1,1);QDCoords(2,1)]; else xi=[QDCoords(1,seg);QDCoords(2,seg)]; xj=[QDCoords(1,seg+1);QDCoords(2,seg+1)]; end % if nod==5||nod==15 % hold on % plot([xi(1),xj(1)],[xi(2),xj(2)]) % end % normal outward unit vector vec=[xj(1)-xi(1);xj(2)-xi(2)]; Ls=sqrt(vec(1)^2+vec(2)^2); nunit=[vec(2);-vec(1)]/Ls; n = [nunit(1) 0 nunit(2);0 nunit(2) nunit(1)]; [gp,wg]=gauleg(ng,-1,1); J=Ls/2; %left boundary segment if xi(1)<tor && xj(1)<tor continue end %down boundary segment if xi(2)<tor && xj(2)<tor continue end %right boundary segment if abs(xi(1)-Lb)<tor && abs(xj(1)-Lb)<tor continue end %up boundary segment if abs(xi(2)-D)<tor && abs(xj(2)-D)<tor for i=1:ng % loop over gauss points - perform numerical integration yg=0.5*(1-gp(i))*xi(2) + 0.5*(1+gp(i))*xj(2); t=-P; fgp = t*wg(i)*J; fg(lin2) = fg(lin2) + fgp; end continue end % Domain internal boundary for i=1:ng % loop over gauss points - perform numerical integration xg=0.5*(1-gp(i))*xi(1)+0.5*(1+gp(i))*xj(1); yg=0.5*(1-gp(i))*xi(2)+0.5*(1+gp(i))*xj(2); gpos=[xg;yg]; v = Domain(gpos,Coords,ds,numnod); % [phi,dphix,dphiy] = shape(gpos,Coords,v,ds); [phi, dphix, dphiy, dphixx, dphiyy, dphixy] = mls_2d(gpos,Coords,v,ds,m); ns = length(v); for j=1:ns col1=(2*v(j))-1; col2=2*v(j); B=[dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)]; kgp = n*Dmat*B*wg(i)*J; kg(lin1,col1) = kg(lin1,col1) + kgp(1,1); kg(lin1,col2) = kg(lin1,col2) + kgp(1,2); kg(lin2,col1) = kg(lin2,col1) + kgp(2,1); kg(lin2,col2) = kg(lin2,col2) + kgp(2,2); end end end end %end loop nodes %IMPOSITION OF K-BOUNDARY CONDITIONS nebc = size(ebc,1); for i=1:nebc nod=ebc(i,1); gpos=[Coords(1,nod);Coords(2,nod)]; v = Domain(gpos,Coords,ds,numnod); [phi, dphix, dphiy, dphixx, dphiyy, dphixy] = mls_2d(gpos,Coords,v,ds,m); ns = length(v); lin1=(nod*2)-1; lin2=nod*2; if ebc(i,2)==1 kg(lin1,:)=0; fg(lin1) = ebc(i,4); for j=1:ns col1=(2*v(j))-1; kg(lin1,col1) = phi(j); end end if ebc(i,3)==1 kg(lin2,:)=0; fg(lin2) = ebc(i,5); for j=1:ns col2=(2*v(j)); kg(lin2,col2) = phi(j); end end end %SOLVE SYSTEM np=kg\fg'; % POST-PROCESSING % Numerical actual displacements [ u ] = getdisp(numnod,Coords,ds,np,m) dlmwrite('ur_patch2.dat',u); % Strain & Stress Computing ps = getstress( numnod,Coords,ds,np,Dmat,m,'rbf_patch'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Plot Nodes for i=1:numnod hold on axis equal str = int2str(i); text(Coords(1,i)+.1,Coords(2,i)-0.1,str) plot(Coords(1,i),Coords(2,i),'o','MarkerFaceColor','k','MarkerSize',3) end