% 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