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