abigNumber=10^12;
%Determine the size of Structural Matrices from Input
Ndofs=Nnodes*Ndim;     %Total Number of Degrees of Freedom


%% Initiallise Structure arrays 
%Stiffness Matrix, Re-arrangment Matrix, Load
%Vector, Displacement Vector

aStructure.Kg=zeros(Ndofs,Ndofs);
aStructure.V=zeros(Ndofs,Ndofs);
aStructure.Pext=zeros(Ndofs,1);
aStructure.U=zeros(Ndofs,1);
aStructure.NBdofs=0;


%% Number the degress of freedom
%Each Structural Node has a set of degrees of freedom that need to be
%atributed with an index
for inode=1:Nnodes
   for idof=1:Ndim
      Nodes{inode}.index(idof)=Ndim*(inode-1)+idof;
   end
end

%% Now Assigned the global degrees of freedom to each one of the elements
%Now assign the Global Dof Numbering to the Elements
for ielement=1:Nelements
   nElementNodes=size(anElement{ielement}.iEnd,2); 
   for inode=1:nElementNodes
    aNode=anElement{ielement}.iEnd(inode); %Node id
        for idof=1:Ndim
            adof=Ndim*(inode-1)+idof;
            anElement{ielement}.index(adof)=Nodes{aNode}.index(idof);  %dof id
        end
    end
end




%% Identify the constrained degrees of freedom
%Get number and iD of consrained dofs
ibdof=0;
for inode=1:Nnodes
   ifix=10^(Ndim-1);
   ifixtemp=Nodes{inode}.ifix;
   if ifixtemp>0
       %then bounded node
       for idof=1:Ndim
           if fix(ifixtemp/ifix)>0
            aStructure.NBdofs=aStructure.NBdofs+1;            
            ifixtemp=ifixtemp-fix(ifixtemp/ifix)*ifix;
            
            %Store Bounded Dof
            adof=Nodes{inode}.index(idof);
            ibdof=ibdof+1;
            aStructure.Bdofs(ibdof,1)=adof;
            
           end
           ifix=fix(ifix/10);
       end
   end   
end