% 
% A program to generate the ground structure for truss topology optimization
%
%Lx=Length of the recatangle (xdir)
%Ly=length of the rectangle (ydir)
%E=Youngs Modu of material used (assuming truss to be made of single material)
%% INPUT Parameters to the program to be given by the user
clear all
clc
Lx = 1;
Ly = 1;
nNodex = 2;
nNodey = 2;
E = 210E9;
%%
Nx=nNodex-1; %Number of nodes along x-axis
Ny=nNodey-1; %Number of nodes along yaxis
                
% Grid node generation
NNODE = 0;
for j = 1:Ny+1,
    for i = 1:Nx+1,
        NNODE = NNODE + 1;
        nx(NNODE) = (Lx/Nx)*(i-1);
        ny(NNODE) = (Ly/Ny)*(j-1);
    end
end
% Grid element generation
NELEM = 0;
for i = 1:NNODE
    for j = i+1:NNODE
        NELEM = NELEM + 1;
        ncon(NELEM,1) = i; 
        ncon(NELEM,2) = j;
    end
end
Ea=ones(NELEM,1)*E;
% Lengths and Young's modulus of elements
for i = 1:NELEM,
    n1 = ncon(i,1);
    n2 = ncon(i,2);
    length(i) = sqrt( (nx(n2)-nx(n1))^2 + (ny(n2)-ny(n1))^2 );

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%Removing Overlapping elements%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ie =1:NELEM,
    eye = ncon(ie,1);
    jay = ncon(ie,2);
    slope(ie) = abs((ny(jay)-ny(eye))/(nx(jay)-nx(eye)));
end
[slpunq,A] = unique(slope);
tol =0.0001;
szslpunq = size(slpunq',1);
for i=1:szslpunq-1,
    if(abs(slpunq(i)-slpunq(i+1))<tol)
        slpunq(i+1) = slpunq(i);
    end
end

slpunq = unique(slpunq);

szslpunq = size(slpunq',1)
distunt = zeros(szslpunq,1);
clc
t = find(slpunq==inf);
tt = find(slope ==inf);
slope(tt)=1e22;
slpunq(t)=1e22;
for i = 1:szslpunq,
    
    for j = 1:NELEM,
        d =length(j);
        if(abs(slope(j)-slpunq(i))<tol)
            
            if abs(length(j)-d)<tol,
                d =length(j);
                distunt(i)=d;
                
            end 
        end
    end
end
ncon(:,4) = zeros(NELEM,1);
p=0;
for i = 1:szslpunq,
    for j = 1:NELEM,
        
        if(abs(slope(j)-slpunq(i))<tol)
            if (abs(length(j) -distunt(i))<tol)
                p=p+1;
                ncon(j,4) = 1;
               
            end
        end
    end
end

r=0;
for(i=1:NELEM)
    if(ncon(i,4)==1)
        r = r+1;
        ncon1(r,1) =ncon(i,1);ncon1(r,2) = ncon(i,2);
    end
end
ncon =ncon1;
NELEM = r;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for(i=1:NELEM)
    Length(i)=sqrt (  (nx(ncon(i,1))-nx(ncon(i,2)))^2 + (ny(ncon(i,1))-ny(ncon(i,2)))^2  );
end

% Plot to test the ground structure

figure(1)
clf
for i = 1:NELEM,
    eye = ncon(i,1); jay = ncon(i,2);
    plot([nx(eye) nx(jay)], [ny(eye) ny(jay)],'-r');
    hold on
end
for(i=1:NNODE)
    plot(nx(i),ny(i),'b*');
    plot(nx(i),ny(i),'bo');
end

for i = 1:NNODE
    xeye = nx(i);
    xjay = ny(i);
    nn = num2str(i);
    text(xeye,xjay,nn,'Fontsize',16);
end
title (['The ground structure'],'Fontsize',20);
hold off
axis off

% Save data to run femtruss.m
save trussdata NNODE NELEM nx ny ncon Length











    


