## spESTtrussLaheWFI.m
close all 
clear
## 
disp('==================================================================')
disp(' spESTtrussLaheWFI.m     ')
disp(' Analysis of a plane truss using the EST method. ')
disp(' The EST method finds the displacements and forces. ')
disp(' The equations matrix are non-symmetric sparse matrix. ')
## http://digi.lib.ttu.ee/opik_eme/slaidid/raamEST1en.pdf
##disp('  ')
disp('---  ')
disp(' http://digi.lib.ttu.ee/opik_eme/slaidid/raamEST1en.pdf ')
disp(' http://digi.lib.ttu.ee/opik_eme/Ehitusmehaanika.pdf#page=408  page=407 ')
## GNU octave version >= 3.0.x
##=========================================================================
## This Program is written by Andres Lahe,   2010-06-26
##                    e-mail: andres.lahe@ttu.ee
## LAST MODIFIED: Andres Lahe,   2013-07-08
## Copyright (c)  2010 by Tallinn University of Technology 
##                Department of Mechanics
##                http://www.ttu.ee/
##
## This program is free software; you can redistribute it and/or
## modify it under the terms of the GNU General Public License
## as published by the Free Software Foundation; either version 2
## of the License, or (at your option) any later version.

## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
## http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
##=========================================================================
##siireVardaA
 eps1=0.000001;
 epsgraf=0.01;
 global spA 
 global siireVardaA
 global SiireJoud
 global X 
# global Ntoerkts
#  global SpTUv
# global SpTWv
disp('---  ')
disp(' daf -- the displacements and the forces at the ends of element. ')
disp(' The element of plane truss have 8 daf''s (4+4). ') 
disp(' The number of unknowns are 8*number_of_elements + number_of_support_reactions. ')
disp(' The matrices are assembled into compressed column sparse matrices. ')
disp('---  ')
disp('  ')
#---------  Number of truss nodes ----------
Number_of_truss_nodes=10
                            SolmedeArv=Number_of_truss_nodes;
#---------  Number of elements ----------
Number_of_elements=17
                            ElementideArv=Number_of_elements;
#---------  Number of supports ----------
Number_of_support_reactions=4 # 4
                            Ntoerkts=Number_of_support_reactions;
disp('  ')

spNNK=8*ElementideArv+Ntoerkts;
   The_number_of_unknowns_are=spNNK
disp('---  ')
#
 d=3.0;  # the panel length
 Npanel=4; #the number of panels
 L=12.0;  # the span of the truss
 H=2.25;  # the height of the truss
 jaotT=L/Npanel; 
#Lp=12.0; % graafika axis 
#
A=9.58E-04  % A=9.58 cm2
%A=1.0
Aa=1.0*A
Ay=1.4*A
Ad=1.2*A
#==========
F=12*1000 # N
Fz=F ;
#==========
disp('--------- Element properties ------------')
#E=2.1E+11  % Elastsusmoodul (E=210GPa) E=2.1E+11 [Pa]
#E=2.1E+05 % kPa
E=2.1E+11 % Pa
#
EAa=Aa*E;
EAy=Ay*E;
EAd=Ad*E;
%EAa =  201180000
%EAy =  2.8165e+08
%EAd =  241416000
disp('--------- Scaling multiplier for displacements----------')
disp(' baasi0=EAy/d; ')   
baasi0=EAy/d;  # scaling multiplier for displacements
#baasi0=EAy/d;
# baasi0=1.0;
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## ---- load variants -----
load_variant=1
#load_variant=2
#load_variant=3
        koormusvariant=load_variant;
#
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch (koormusvariant)
case{1}
disp(' Load variant 1 ')
#
#disp(' Node forces in global  coordinates ')
# sSolmF(forces,1,nodes); forces=[Fx; Fz]
sSolmF = zeros(2,1,SolmedeArv);
#sSolmF(:,1,1)=0.0
sSolmF(2,1,2)=Fz/2;  #  s2F(2,1)=F/2;
#sSolmF(:,1,3)=0.0
sSolmF(2,1,4)=Fz;    #   s4F(2,1)=F;
#sSolmF(:,1,5)=0.0
sSolmF(2,1,6)=Fz;    #   s6F(2,1)=F;
#sSolmF(:,1,7)=0.0
sSolmF(2,1,8)= Fz;   #   s8F(2,1)=F;
#sSolmF(:,1,9)=0.0
sSolmF(2,1,10)=Fz/2;  #  s10F(2,1)=F/2;
#
# Support shift - tSiire#
# Support shift multiplied by scaling multiplier
tSiire = zeros(2,1,SolmedeArv);
#tSiire(:,1,1)= 0.0
#tSiire(2,1,1)= 0.01*baasi0
# Support shift - tSiire#
# Support shift multiplied by scaling multiplier
tSiire = zeros(3,1,SolmedeArv);
#tSiire(:,1,1)= 0.0
#tSiire(2,1,1)= 0.01*baasi0
#tSiire;
##%%%%
case{2}
disp(' Load variant 2 ')
case{3}
disp(' Load variant 3 ')
#%%%%%
otherwise
disp(' No load variant cases ')
endswitch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('------  ')
disp(' Node forces in global  coordinates ')
disp(' Fx Fz; 1 ; Node number; ')
disp('------  ')
sSolmF(:,1,:)=sSolmF(:,1,:)
#
#==========
#     Nodal coordinates
#==========
krdn=[#	 x       z
        0.0      0.0  ;   # node 1
        0.0     -2.25 ;   # node 2
        3.0     -0.75 ;   # node 3
        3.0     -2.25 ;   # node 4
        6.0     -1.5  ;   # node 5
        6.0     -2.25 ;   # node 6
        9.0     -0.75 ;   # node 7
        9.0     -2.25 ;   # node 8
       12.0      0.0  ;   # node 9
       12.0     -2.25];   # node 10
#========== 
#
#==========
#  Restrictions on support displacements.
# Support  No   u   w  fi   hold on - 1, open - 0
#==========       
tsolm=[1    1 1 0;   # node 1
       9    1 1 0];  # node 9
%==========
TSARV=size(tsolm);
NTSARV=TSARV(1,1); 
##%%%%%%%%%%%
ToeSidemeteArv=0;
for i=1:NTSARV
  for j=2:3
     if (tsolm(i,j) == 1)
       ToeSidemeteArv=ToeSidemeteArv+1;
     endif 
  endfor
endfor
#
#%%%
#ToeSidemeteArv
if ((Ntoerkts-ToeSidemeteArv)  !=0 )
disp(' ')
disp(' ------- ')
disp(' ')
Number_of_support_reactions
Error_by_inserting_support_reactions=sprintf('Number of support reactions must be, and is  %d, %d',Ntoerkts,ToeSidemeteArv)
disp(' ')
error(' Error inserting the restrtictions on the support displacements.')
endif
#
#%%%
#==========
# ------------- Element topology, properties  ---------
#==========
elasts=[# Element properties, 
###selemjl=[# Element properties,
#        n2 - end of the element, 
#            n1 - beginning of the element, 
#                 A  E;  
          2   1  Ad  E;  # element 1
          3   1  Aa  E;  # element 2
          3   2  Ad  E;  # element 3
          4   2  Ay  E;  # element 4
          4   3  Ad  E;  # element 5
          6   4  Ay  E;  # element 6
          5   4  Ad  E;  # element 7
          5   3  Aa  E;  # element 8
          6   5  Ad  E;  # element 9
          8   6  Ay  E;  # element 10
          8   5  Ad  E;  # element 11
          7   5  Aa  E;  # element 12
          8   7  Ad  E;  # element 13
         10   8  Ay  E;  # element 14
         10   7  Ad  E;  # element 15
          9   7  Aa  E;  # element 16
         10   9  Ad  E]; # element 17
#
#==========
%%selemjl(:,4)=E; % E -elastsusmoodul
%
%selemE=[selemjl(:,1) selemjl(:,2) selemjl(:,3) selemjl(:,4)];
EARV=size(elasts);
NEARV=EARV(1,1);

%disp('--------------------------------------------------------------')
%
%==========
SARV=size(krdn);
NSARV=SARV(1,1);
%%%%%%%%%%%%%
for i=1:NSARV
sJrN(i)=i;
end
%
sJrNT=sJrN';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%disp('-----------------------------')
disp('===============================')
disp('    Nodal coordinates ')
disp('  Node     X        Z   ')
disp('-----------------------------')
%disp([krdn])
for i=1:NSARV
disp(sprintf('  %2i    %7.4f  %7.4f',sJrNT(i),    krdn(i,1), krdn(i,2)))
endfor
disp('-----------------------------')
%%%%
if ((SolmedeArv-NSARV) != 0)
disp(' ')
disp(' ------- ')
disp(' ')
Number_of_truss_nodes
Error_by_nodes_input=sprintf('Number of nodes must be, but is: %d, %d', SolmedeArv, NSARV)
disp(' ')
error(' Error by nodes input.')
endif
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========
TSARV=size(tsolm);
NTSARV=TSARV(1,1);
%%%%%%%%%%%%%
for i=1:NTSARV
sJrN(i)=i;
end
%
sTJrNT=sJrN';
%%%%%%%%%%%%
%tsolm
%disp('--------------------------------')
disp('==================================================')
disp(' Restrictions on support displacements. ')
disp(' No    Node No    u   w  ;  hold on-1, open-0')
disp('--------------------------------------------------')
%disp([tsolm])
for i=1:NTSARV
disp(sprintf('  %2i     %2i      %2i  %2i  %2i ',sTJrNT(i),    tsolm(i,1), tsolm(i,2), tsolm(i,3), tsolm(i,4)))
endfor
disp('---------------------------------')
%%%
TOEARV=size(tsolm);
TOARV=TOEARV(1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#%%
# ------------- Element properties, topology and hinges ---------
#-----------------------
slmnr=0;
for i=1:NEARV
  for j=1:8
   slmnr=slmnr+1;
   selemjln(i,j)=slmnr;
  endfor
endfor
#
selemjl=[selemjln elasts];
#disp('===============================================================================================')
#disp('  Displacements and the forces numbers (DaFs) ')
#disp('  at the end and at the beginning of elements.          Nodes                  ')
#disp('       At the end:              At the beginning:                                              ')
#disp('    u      w     fi    N     u     w    fi    N     The_end The_beginning        A       E     ')
#disp('-----------------------------------------------------------------------------------------------')
##selemjln;
##
#for i=1:NEARV
#disp(sprintf('  %3i    %3i   %3i   %3i   %3i   %3i   %3i   %3i      %3i     %3i          %9.3e   %9.3e', #selemjl(i,1), ...
#selemjl(i,2), selemjl(i,3),selemjl(i,4), ...
#selemjl(i,5), selemjl(i,6), selemjl(i,7),selemjl(i,8), ...
#selemjl(i,9), selemjl(i,10), selemjl(i,11),selemjl(i,12)))
#endfor

Displacements_and_the_forces_numbers=selemjl;
#
#selemjl=[selemjln elasts];
selemE=[selemjl(:,1:12)];
#EARV=size(selemE);
#NEARV=EARV(1,1);Jr
#%%%%%%%%%%%%
for i=1:NEARV
JrN(i,1)=i;
endfor
#
JrNT=JrN';
#%%%%%%%%%%
disp('===============================================================================================')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%selemE=[selemjl(:,1) selemjl(:,2) selemjl(:,3) selemjl(:,4)];
EARV=size(selemjl);
NEARV=EARV(1,1);
for i=1:NEARV
selemEA(i,1)=elasts(i,3)*elasts(i,4);
endfor
%EA=selemEA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:NEARV
 selem(i,1:11) = [selemE(i,1:10)  selemEA(i,1)];
endfor
#
#selem
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%disp(' olen real 281   ')
%pause
%for i=1:NEARV+3
for i=1:NEARV
JrN(i,1)=i;
end
%17
JrNT=JrN';
%%%%%%%%%%%%
#%disp('-----------------------------------')
#disp('====================================')
#disp(' Axial stiffness of element ')
#disp('   No        EA          ')
#disp('------------------------------------')
#%disp([krdn])
#for i=1:NEARV
#disp(sprintf('  %2i    %4.3e',JrN(i),    selemEA(i,1)))
#endfor
#disp('--------------------------------------')
%%%%%%
%for i=1:17
%  for j=1:8
%   n=n+1;
%   Ns(i,j)=n;
%  endfor
%endfor
#%srtkElem=[
#     1     2     3     4     5     6     7     8;
#     9    10    11    12    13    14    15    16;
#    17    18    19    20    21    22    23    24;
#    25    26    27    28    29    30    31    32;
#    33    34    35    36    37    38    39    40;
#    41    42    43    44    45    46    47    48;
#    49    50    51    52    53    54    55    56;
#    57    58    59    60    61    62    63    64;
#    65    66    67    68    69    70    71    72;
#    73    74    75    76    77    78    79    80;
#    81    82    83    84    85    86    87    88;
#    89    90    91    92    93    94    95    96;
#    97    98    99   100   101   102   103   104;
#   105   106   107   108   109   110   111   112;
#   113   114   115   116   117   118   119   120;
#   121   122   123   124   125   126   127   128;
#   129   130   131   132   133   134   135   136];
%%%
%selem=[selemjl(:,1) selemjl(:,2) selemjl(:,3) selemjl(:,4)];
%selem=[selemjl(:,1) selemjl(:,2)     srtkElem(:,1:8)];
selem;
%%%%%%%%%%%%%
#for i=1:NEARV+3
for i=1:NEARV
JrN(i,1)=i;
end
%
JrNT=JrN';
%%%%%%%%%%%%
disp('===============================================================================================')
disp('                           Topology of elements ')
disp('                                   At the end:               At the beginning:                 ')
disp(' No       end   beginning      u     w     fi    N        u     w     fi      N          EA    ')
disp('-----------------------------------------------------------------------------------------------')
%disp([selem])
for i=1:NEARV
disp(sprintf('  %2i     %3i   %3i           %3i   %3i    %3i  %3i       %3i   %3i    %3i    %3i      %4.3e',JrN(i,1), selem(i,9), selem(i,10),   selem(i,1), selem(i,2), selem(i,3),selem(i,4), ...
selem(i,5), selem(i,6), selem(i,7),selem(i,8), ...
 selem(i,11)))
endfor
disp('------------------------------------------------------------------------------------------')
%%%%%%%%%%%%
EARV=size(selem);
NEARV=EARV(1,1);
#%%%
if ((ElementideArv-NEARV)  !=0 )
disp(' ')
disp(' ------- ')
disp(' ')
Number_of_elements
Error_by_elements_input=sprintf(' Number of elements must be, and is  %d, %d',ElementideArv,NEARV)
disp(' ')
error(' Error inserting the topology of elements.')
endif
#%%%
#% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#for i=1:NSARV
#LkoordN=selem(i,9);
#AkoordN=selem(i,10);
#Koordinaadid=[LkoordN AkoordN]
#endfor

#selem
lvarras=VardaPikkusTr(NSARV,NEARV,krdn,selem);
#lvarras
#%%%%%%%%%%%
disp('   ')
NEARV;
NNK=8*NEARV+Ntoerkts;

lvarras;
#
#selemjl
#%disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ')
for i=1:NEARV
elemendiN(i,1)=i;
siireVardaA(i,1:3)=selemjl(i,5:7);
siireVardaL(i,1:3)=selemjl(i,1:3);
joudVardaA(i,1)=selemjl(i,8);
joudVardaL(i,1)=selemjl(i,4);
endfor
#
siireVardaA;
siireVardaL;
joudVardaA;
joudVardaL;
#%%%%%%%%%%%%

%%%%%ABB=ABB
%NNK;

#--------- Define a matrix spA to be sparse -----------
disp('-------- Sparse matrix instantiation --------  ')
disp('  spA=sparse(NNK,NNK)  ')
spA=sparse(NNK,NNK) % nullistan võrrandisüsteemi kordajate maatriksi
#--------- Define a vector B ------------
disp('-------- Right-hand side of the equations (RHS). --------  ')
disp('  B=zeros(NNK,1);  ')
B=zeros(NNK,1);   % nullistan võrrandisüsteemi vabaliikmed
disp('   ')
%%%%%%%%%%%%%%%%%%%%%%%%%



#------------ The basic equations of frame ----------------
#--------- Define a matrix spA to be sparse -----------
disp('  ')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ')
disp(' function AlgParm=LaheTrussDFI(baasi0,Ntoerkts,sSolmF,tsolm,tSiire,krdn,selem) ')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ')
disp('  ')
#AlgParm=LaheTrussDFIpoolik(baasi0,Ntoerkts,sSolmF,tsolm,tSiire,krdn,selem);
AlgParm=LaheTrussDFI(baasi0,Ntoerkts,sSolmF,tsolm,tSiire,krdn,selem);
disp('  ')
AlgPar=AlgParm;
disp('  ')
%
#disp('=============================================================')
#disp('  Unscaled initial parmeter vector (Sign Convention 2)   ')
#disp('Element No    u          w          fi             N      ')
#disp('--------------------------------------------------------------')
#for i=1:NEARV
#disp(sprintf('  %3i     %9.3e   %9.3e   %9.3e    %9.3f',JrN(i,1),  AlgPar(i,1),   AlgPar(i,2),   #AlgPar(i,3),  AlgPar(i,4)))
##%AlgPar
#endfor

#disp('---------------------------------------------------------------')
%
NEARV2=2*NEARV;
ElemNr=0;
for i=1:NEARV
  for j=1:2
      ElemNr=ElemNr+1;
      NrElem(ElemNr,1)=i;
   endfor
endfor
%NrElem'
%
disp('============================================================================')
disp('   Unscaled displacements and forces at beginning / at end  ')
disp('                   (local coordinates, Sign Convention 1) ')
disp('  No  Element_No      u          w          fi             N  ')
disp('----------------------------------------------------------------------------')
for i=1:NEARV2
JrN2(i,1)=i;
disp(sprintf('  %3i   %3i     %9.3e   %9.3e   %9.3e    %11.3f',JrN2(i,1), NrElem(i,1), SiireJoud(i,1),   SiireJoud(i,2),   SiireJoud(i,3),  SiireJoud(i,4)))
endfor
disp('----------------------------------------------------------------------------')
disp('============================================================================')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:NEARV
LkoordN=selem(i,9);
AkoordN=selem(i,10);
#LkoordN=selemjl(i,1);
#AkoordN=selemjl(i,2);
DeltaX(i)=krdn(LkoordN,1)-krdn(AkoordN,1);
DeltaZ(i)=krdn(LkoordN,2)-krdn(AkoordN,2);
VGRx(i,1)=krdn(AkoordN,1);
VGRx(i,2)=krdn(LkoordN,1);
VGRz(i,1)=krdn(AkoordN,2);
VGRz(i,2)=krdn(LkoordN,2);
endfor
%
%suurused=['      siire  u - ';
%          '      siire  w - ';
%          '      pööre fi - ';
%          ' Normaaljõud N - ';
%%%%%%%%%%%%%
figure(1)
ax = gca ();
axy=[ 0   0   1   0.55];
set (ax, "outerposition", [ 0   0   1   0.55])
%
grid off
%
Lkorg=L;
%
for j=1:NEARV
       IR=j;
 % 
   axis([-d L+d -H-2.0 1.0], "ij")
 %
plot(VGRx(IR,:),VGRz(IR,:),"3")
hold on
endfor
%
%a1=eval(char([194 176]))
%title('S(char([194 176]))restikskeem',"fontsize",12)
title('spESTtrussLaheWFI',"fontsize",12)
ax=gca();
      xlabel('x',"fontsize",12)
      ylabel('z',"fontsize",12)
 xticks(1,:)=0.0:jaotT:L;
 set (ax, "xtick", xticks)
%
%%Sõlmede numbrid
TTT=' ';
TST=' ';
TT1=' ';
TT1=' ';
%%%%%%%%%%
for i=1:NSARV;
%
   IR=i;
   TTT=eval(sprintf('%i',IR));
   TST=int2str(TTT);
 
text(krdn(IR,1), krdn(IR,2)+0.17, TST)
#
endfor
#
# Varraste numbrid
for j=1:NEARV;
#
   IR=j;
   TTT=eval(sprintf('%i',IR));
   TST=int2str(TTT);
text((VGRx(IR,1)+VGRx(IR,2))/2+0.1, (VGRz(IR,1)+VGRz(IR,2))/2, TST)
#
endfor
#text(L/2-d/2, -H-d/2, 'spESTtrussLaheWFI')
#
hold off
#%%%%%%
#axis([-d L+d 1.0 -H-2.0], "ij")
#pause (1)
#
#%Vaata: http://www.gnu.org/software/octave/doc/interpreter/Printing-Plots.html
#
#print('-deps','spESTtrussLaheWFI.eps');
print('-dpng','spESTtrussLaheWFI.png');
#print('-dfig','-landscape','-mono','-solid','spESTtrussLaheWFI.fig');
refresh
## %%%%%%%%%%%%%%%%%%%%%
disp(' Support reactions ')
C1=X(137,1)
C2=X(138,1)
C3=X(139,1)
C4=X(140,1)
disp('---------------')
##
##
N5=X(36,1)
H1=X(137,1)
Va=X(138,1)

disp('---------------------------------')
## N5=-13.71305
## H=((24-6)*9+N5*6)/2.25
##disp(sprintf('%11.6e     %11.6e      %11.6e ', N5,   H1,   Va ))
disp(' Testing a static equilibrium for the truss. ')
disp(' Sum of the moments of forces acting about point 8, M8. ')
disp(' The elements  4, 5, 8 are cutted through.  ')  
##disp('  Momentpunktiks on valitub sõlm 8   ') 
disp('  ')
disp('  M8=(Va+6000.0)*9.0-N5*6.0+H1*2.25  ') 
 M8=(Va+6000.0)*9.0-N5*6.0+H1*2.25 
disp('  ')
disp('-----------------------------------')
#%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
spy(spA)
title('spy(spA,14) - the sparse matrix spA(140,140) non zero elements [1.9%] ')
refresh
#
#print('-deps','pESTtruss_sparse_matrix.eps');
#print('-dpng','pESTtruss_sparse_matrix.png');
print('-dfig','-landscape','spESTtruss_sparse_matrix.fig');
##print('-dfig','-landscape','-mono','-solid','spESTtruss_sparse_matrix');
refresh
#
## Compressed Column Sparse (rows = 140, cols = 140, nnz = 376 [1.9%])
#
# octave
#disp('-----------------------------------')
disp('  ')
disp(' The END ')
#%%%%%%%%%%%%%%%

