## spESTGerberBeamWFI.m
close all 
clear
## 
disp('==================================================================')
disp(' spESTGerberBeamWFI.m     ')
disp(' Analysis of a Gerber beam using the EST method. ')
disp(' The EST method finds the displacements and forces. ')
disp(' The equations matrix are non-symmetric sparse matrix. ')
disp('  ')
disp(' http://digi.lib.ttu.ee/opik_eme/Ehitusmehaanika.pdf#page=437  #page=269           ')
## GNU octave version >= 3.0.x
##=========================================================================
## This Program is written by Andres Lahe,   2010-07-16
##                    e-mail: andres.lahe@ttu.ee
## LAST MODIFIED: Andres Lahe,   2010-08-03
## 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
##=========================================================================
##
 eps1=0.000001;
 epsgraf=0.01;
 global spA
 global X  
 global siireVardaA
disp('---  ')
disp(' daf -- the displacements and the forces at the ends of element. ')
disp(' The element of a continious beam 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 frame nodes ----------
Number_of_beam_nodes=8  #5
                            SolmedeArv=Number_of_beam_nodes;
#---------  Number of elements ----------
Number_of_elements=7    #4
                            ElementideArv=Number_of_elements;
#---------  Number of supports ----------
Number_of_support_reactions=5
                            Ntoerkts=Number_of_support_reactions;
disp('  ')

NNK=8*ElementideArv+Ntoerkts;
   The_number_of_unknowns_are=NNK
disp('---  ')
disp(' The displacements and the forces calculated on the parts of the element ''Nmitmeks''.  ')
#Nmitmeks=4
Nmitmeks=5
d=10;
#Lp=Ls; # graphics axis
#==========
disp('--------- Element properties ------------') 
#
#E=2.1E+11  # Elastsusmoodul (E=210GPa) E=2.1E+11 [Pa]
#E=2.1E+05 # kPa
## I-40  Hristloige=400mm=0.4m   I=1902cm4=1.902E+03cm4=1.902E-05m4
#
EI=20000
GA=1.0*10^15;
%F1=60;
%F2=40;
%F3=80;
%qz=12.0;
#qz=0.0;
%
F1s=20.0;
F2k=40.0;
aF2k=4.0;
F3k=30.0;
aF3k=2.0;
F4k=100.0;
aF4k=4.0;
F5s=20.0;
%
qk=8;
%
L1k=4.0
L2k=8.0
L3k=2.0
L4k=6.0
L5k=2.0
L6k=8.0
L7k=1.0;
Ls=(L1k+L2k+L3k+L4k+L5k+L6k+L7k);
Lp=Ls;
baasi0=EI/8.0
##baasi0=1.0; 
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## load variants
load_variant=1
#load_variant=2
#load_variant=3
        koormusvariant=load_variant
#koormusvariant=2
#koormusvariant=3
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch (koormusvariant)
case{1}
#
disp(' The load variant 1 ')
#disp(' Element load  in local coordinates ')
# The uniformly distributed load in local coordinate z direction
LoadsqONelement=4;
esQkoormu=zeros(LoadsqONelement,3,ElementideArv);
esQkoormus(1,1:3,1)=[0.0 0.0 L1k];
esQkoormus(1,1:3,2)=[qk 0.0 L2k];
esQkoormus(1,1:3,3)=[qk 0.0 L3k];
esQkoormus(1,1:3,4)=[qk 0.0 L4k];
esQkoormus(1,1:3,5)=[0.0 0.0 L5k];
esQkoormus(1,1:3,6)=[0.0 0.0 L6k];
esQkoormus(1,1:3,7)=[0.0 0.0 L7k];
esQkoormus
#
# The point load in local coordinate z direction
LoadsF_on_Element=5; 
# 2 - Fz, aFz - the coordinate of the point of force Fz application
esFjoud=zeros(LoadsF_on_Element,2,ElementideArv);
esFjoud(1,1:2,1)=[0.0  L1k];
esFjoud(1,1:2,2)=[F2k aF2k];  
esFjoud(1,1:2,3)=[0.0  L3k];
esFjoud(1,1:2,4)=[F3k aF3k];
esFjoud(1,1:2,5)=[0.0  L5k];
esFjoud(1,1:2,6)=[F4k aF4k]; 
esFjoud(1,1:2,7)=[0.0  L7k];
#esFjoud
#disp(' Node forces in global  coordinates ')
#disp(' Node 1 load ')
# sSolmF(forces,1,nodes); forces=[Fz; My]
sSolmF = zeros(2,1,SolmedeArv);
#sSolmF(:,1,1)= 0.0
sSolmF(1,1,2)= F1s;  # Fz
sSolmF(2,1,2)= 0.0;  # My
#sSolmF(:,1,3)= 0.0
#sSolmF(:,1,4)= 0.0
#sSolmF(:,1,5)= 0.0
#sSolmF(:,1,6)= 0.0
#sSolmF(:,1,7)= 0.0
sSolmF(1,1,8)= F5s;   # Fz
sSolmF(2,1,8)= 0.0;   # My
%disp(' Kustuta ')
     s1F(:,:) = sSolmF(:,1,1);
     s2F(:,:) = sSolmF(:,1,2);
     s3F(:,:) = sSolmF(:,1,3);
     s4F(:,:) = sSolmF(:,1,4);
     s5F(:,:) = sSolmF(:,1,5);
     s6F(:,:) = sSolmF(:,1,6); 
     s7F(:,:) = sSolmF(:,1,7);
     s8F(:,:) = sSolmF(:,1,8);
#
#s1F(1,1)=0.0;  # force Fz
#s1F(2,1)=0.0;    # force My
#%%%%%
# The support shift - tSiire
# The support shift multiplied by scaling multiplier
tSiire = zeros(2,1,SolmedeArv);
#tSiire(:,1,1)= 0.0
#tSiire(2,1,1)= 0.01*baasi0
#tSiire(:,1,2)= 0.0
#tSiire(:,1,3)= 0.0
#tSiire(:,1,4)= 0.0
#tSiire(:,1,5)= 0.0
##%%%%
case{2}
disp(' The load variant 2 ')
#disp(' Element load  in local coordinates ')
# The uniformly distributed load in local coordinate z direction
LoadsqONelement=4;
esQkoormu=zeros(LoadsqONelement,3,ElementideArv);
esQkoormus(1,1:3,1)=[0.0 0.0 0.0];
esQkoormus(1,1:3,2)=[0.0 0.0 0.0];
esQkoormus(1,1:3,3)=[0.0 0.0 0.0];
esQkoormus(1,1:3,4)=[0.0 0.0 0.0];
esQkoormus(1,1:3,5)=[0.0 0.0 0.0];
esQkoormus(1,1:3,6)=[0.0 0.0 0.0];
esQkoormus(1,1:3,7)=[0.0 0.0 0.0];
esQkoormus
# The point load in local coordinate z direction
LoadsF_on_Element=5; 
# 2 - Fz, aFz - the coordinate of the point of force Fz application
esFjoud=zeros(LoadsF_on_Element,2,ElementideArv);
esFjoud(1,1:2,1)=[0.0 0.0];
esFjoud(1,1:2,2)=[0.0 0.0];
esFjoud(1,1:2,3)=[0.0 0.0];  
esFjoud(1,1:2,4)=[0.0 0.0];
esFjoud(1,1:2,5)=[0.0 0.0];
esFjoud(1,1:2,6)=[0.0 0.0];
esFjoud(1,1:2,7)=[0.0 0.0];
esFjoud
#
#disp(' Node forces in global  coordinates ')
#disp(' Node 1 load ')
# sSolmF(forces,1,nodes); forces=[Fz; My]
sSolmF = zeros(2,1,SolmedeArv);
#sSolmF(:,1,1)= 0.0
#sSolmF(:,1,2)= 0.0
#sSolmF(:,1,3)= 0.0
#sSolmF(:,1,4)= 0.0
#sSolmF(:,1,5)= 0.0
#sSolmF(:,1,6)= 0.0
#sSolmF(:,1,7)= 0.0
sSolmF(1,1,8)= F5s;   # Fz
sSolmF(2,1,8)= 0.0;   # My
%disp(' Kustuta ')
     s1F(:,:) = sSolmF(:,1,1);
     s2F(:,:) = sSolmF(:,1,2);
     s3F(:,:) = sSolmF(:,1,3);
     s4F(:,:) = sSolmF(:,1,4);
     s5F(:,:) = sSolmF(:,1,5);
     s6F(:,:) = sSolmF(:,1,6); 
     s7F(:,:) = sSolmF(:,1,7);
     s8F(:,:) = sSolmF(:,1,8);
#
#s1F(1,1)=0.0;  # force Fz
#s1F(2,1)=0.0;    # force My
#%%%%%
# The support shift - tSiire
# The support shift multiplied by scaling multiplier
tSiire = zeros(2,1,SolmedeArv);
#tSiire(:,1,1)= 0.0
#tSiire(2,1,1)= 0.01*baasi0
#tSiire(:,1,2)= 0.0
#tSiire(:,1,3)= 0.0
#tSiire(:,1,4)= 0.0
#tSiire(:,1,5)= 0.0
##%%%%
case{3}
disp(' The load variant 3 ')
#disp(' Element load  in local coordinates ')
# The uniformly distributed load in local coordinate z direction
LoadsqONelement=4;
esQkoormu=zeros(LoadsqONelement,3,ElementideArv);
esQkoormus(1,1:3,1)=[0.0 0.0 0.0];
esQkoormus(1,1:3,2)=[0.0 0.0 0.0];
esQkoormus(1,1:3,3)=[0.0 0.0 0.0];
esQkoormus(1,1:3,4)=[0.0 0.0 0.0];
esQkoormus(1,1:3,5)=[0.0 0.0 0.0];
esQkoormus(1,1:3,6)=[0.0 0.0 0.0];
esQkoormus(1,1:3,7)=[0.0 0.0 0.0];
esQkoormus
# The point load in local coordinate z direction
LoadsF_on_Element=5; 
# 2 - Fz, aFz - the coordinate of the point of force Fz application
esFjoud=zeros(LoadsF_on_Element,2,ElementideArv);
esFjoud(1,1:2,1)=[0.0 0.0];
esFjoud(1,1:2,2)=[0.0 0.0];
esFjoud(1,1:2,3)=[0.0 0.0];  
esFjoud(1,1:2,4)=[0.0 0.0];
esFjoud(1,1:2,5)=[0.0 0.0];
esFjoud(1,1:2,6)=[0.0 0.0];
esFjoud(1,1:2,7)=[0.0 0.0];
esFjoud
#
#disp(' Node forces in global  coordinates ')
#disp(' Node 1 load ')
# sSolmF(forces,1,nodes); forces=[Fz; My]
sSolmF = zeros(2,1,SolmedeArv);
#sSolmF(:,1,1)= 0.0
#sSolmF(:,1,2)= 0.0
#sSolmF(:,1,3)= 0.0
#sSolmF(:,1,4)= 0.0
#sSolmF(:,1,5)= 0.0
#sSolmF(:,1,6)= 0.0
#sSolmF(:,1,7)= 0.0
sSolmF(1,1,8)= F5s;   # Fz
sSolmF(2,1,8)= 0.0;   # My
%disp(' Kustuta ')
     s1F(:,:) = sSolmF(:,1,1);
     s2F(:,:) = sSolmF(:,1,2);
     s3F(:,:) = sSolmF(:,1,3);
     s4F(:,:) = sSolmF(:,1,4);
     s5F(:,:) = sSolmF(:,1,5);
     s6F(:,:) = sSolmF(:,1,6); 
     s7F(:,:) = sSolmF(:,1,7);
     s8F(:,:) = sSolmF(:,1,8);
#
#s1F(1,1)=0.0;  # force Fz
#s1F(2,1)=0.0;    # force My
#%%%%%
# The support shift - tSiire
# The support shift multiplied by scaling multiplier
tSiire = zeros(2,1,SolmedeArv);
#tSiire(:,1,1)= 0.0
#tSiire(2,1,1)= 0.01*baasi0
#tSiire(:,1,2)= 0.0
#tSiire(:,1,3)= 0.0
#tSiire(:,1,4)= 0.0
#tSiire(:,1,5)= 0.0
##%%%%
otherwise
disp(' No load variant cases ')
endswitch
%%%%%%%%%%%%%%
disp(' The element loads  in local coordinates ')
#qz=qzZ'
#aL=aLXx'
#Fz=FZz'
esQkoormus
esFjoud

disp(' Node 1 load in global coordinates ')
s1F'
disp(' Node 2 load in global coordinates ')
s2F'
disp(' Node 3 load in global coordinates ')
s3F'
disp(' Node 4 load in global coordinates ')
s4F'
disp(' Node 5 load in global coordinates ')
s5F'
disp(' Node 5 load in global coordinates ')
s6F'
disp(' Node 5 load in global coordinates ')
s7F'
disp(' Node 5 load in global coordinates ')
s8F'
#
#==========
# Nodal coordinates
#==========
krdn=[	0.0  ;	# node 1
	4.0  ;	# node 2
       12.0  ;	# node 3
       14.0  ;	# node 4
       20.0  ;	# node 5
       22.0  ;	# node 6
       30.0  ;	# node 7
       31.0 ];	# node 8   
#==========
#
#==========
#  The restrtictions on the support displacements.
# Node Nr   w  fi   hold on - 1, open - 0
#==========
tsolm=[1    1 1  ;	# node 1
       3    1 0  ;	# node 3
       6    1 0  ;	# node 6
       7    1 0  ];	# node 7
#==========
# The support shift - tSiire is enterd in load variants
# The support shift multiplied by scaling multiplier
#==========
# ------------- Element properties, topology and hinges ---------
#==========
elasts=[# Element properties, 
#                     n2 - end of the element, 
#                          n1 - beginning of the element,
#                               Q, M hinges at end of the element
#                                        Q, M hinges at beginning of the element. 
        EI  GA    2     1    0 1    0 0;    % element 1
	EI  GA    3     2    0 0    0 1;    % element 2
	EI  GA    4     3    0 1    0 0;    % element 3
	EI  GA    5     4    0 1    0 1;    % element 4
	EI  GA    6     5    0 0    0 1;    % element 5
	EI  GA    7     6    0 0    0 0;    % element 6	
	EI  GA    8     7    0 0    0 0];   % element 7	
# 1-hinge 'true' (shear-, moment hinge)
#==========
#
# Validations
#
##%%%%%%%%%%%
SARV=size(krdn);
NSARV=SARV(1,1);
##%%%%%%%%%%%
for i=1:NSARV
sJrN(i)=i;
end
#
sJrNT=sJrN';
##%%%%%%%%%%
#disp('-----------------------------')
disp('===============================')
disp(' Nodal coordinate ')
disp('  Node    X           ')
disp('-----------------------------')
#disp([krdn])
for i=1:NSARV
disp(sprintf('  %2i    %7.4f  %7.4f',sJrNT(i),    krdn(i,1)))
endfor
disp('-----------------------------')
##%%
if ((SolmedeArv-NSARV) != 0)
disp(' ')
disp(' ------- ')
disp(' ')
Number_of_beam_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);
##%%%%%%%%%%%
ToeSidemeteArv=0;
for i=1:NTSARV
sJrN(i)=i;
  for j=2:3
     if (tsolm(i,j) == 1)
       ToeSidemeteArv=ToeSidemeteArv+1;
     endif 
  endfor
endfor
#
sTJrNT=sJrN';
##%%%%%%%%%%
%ToeSidemeteArv;
#disp('--------------------------------')
disp('================================================')
disp(' The restrtictions on the support displacements.')
disp('JrkNr  Node Nr  w  fi;  hold on-1, open-0 ')
disp('------------------------------------------------')
#disp([tsolm])
for i=1:NTSARV
disp(sprintf('  %2i     %2i    %2i %2i ',sTJrNT(i),    tsolm(i,1), tsolm(i,2), tsolm(i,3)))
endfor
disp('------------------------------------------------')
##%%
##%%%%%%%%%%%
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 properties, topology and hinges ---------
#
disp('===============================================')
disp('        Topology and  hinges. ')
disp('  End of the element, beginning of the element,  ')
%disp('  shear-, moment hinge; 1 - hinge ''true'' ')
disp('-----------------------------------------------')
elemkinni=elasts(:,3:8);
Topology_and_and_hinges=elemkinni
EARV=size(elemkinni);
NEARV=EARV(1,1);
NEARV;
#%%%
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
#%%%
#%%%%%%%%%%%%%%%%%%%%
# --------The displacements and the forces numbers at the end and beginning of element.  ----- -----
slmnr=0;
for i=1:NEARV
 for j=1:8
   slmnr=slmnr+1;
   selemjln(i,j)=slmnr;
 endfor
endfor
#
disp('====================================================')
disp('  The displacements and the forces numbers (daf''s) ')
disp('  at the end and at the beginning of elements.  ')
disp('----------------------------------------------------')
selemjln;
The_displacements_and_the_forces_numbers=selemjln
#
selemjl=[selemjln elasts];
#==========
#selemjl=[1  2  3  4  5  6  7  8 EI GA 2 1 ;   # element 1
#         9 10 11 12 13 14 15 16 EI GA 3 2 ;   # element 2
#	 17 18 19 20 21 22 23 24 EI GA 4 3 ;   # element 3
#	 25 26 27 28 29 30 31 32 EI GA 5 4 ];  # element 4
#==========
selem=[selemjl(:,1:16)];
#%%%%%%%%%%%%
lvarras=VardaPikkusT(NSARV,NEARV,krdn,selem);
#%%%%%%%%%%%
for i=1:NEARV
JrN(i,1)=i;
endfor
#
JrNT=JrN';
#%%%%%%%%%%%
disp('=======================================================================================================')
disp(' Topology of elements                                                                      Moment hinge')
disp('        At the end:             At the beginning:                            Node number   ''true''=1  ')
disp(' Nr.    w,   fi,   Q,   M.     w,   fi,   Q,   M.       EI            GA     End Beginning   M    M')
disp('-------------------------------------------------------------------------------------------------------')
#
for i=1:NEARV
disp(sprintf(' %2i     %2i   %2i   %2i   %2i      %2i   %2i   %2i   %2i     %9.3e    %9.3e  %2i  %3i      %4i   %2i',JrN(i,1),    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), ...
selemjl(i,13), selemjl(i,14)))
endfor
disp('-------------------------------------------------------------------------------------------------------')
#%%%
EARV=size(selem);
NEARV=EARV(1,1);
%%selemjl
for i=1:NEARV
elemendiN(i,1)=i;
siireVardaA(i,1:2)=selemjl(i,5:6);
siireVardaL(i,1:2)=selemjl(i,1:2);
joudVardaA(i,1:2)=selemjl(i,7:8);
joudVardaL(i,1:2)=selemjl(i,3:4);
endfor
#
#--------- Define a matrix spA to be sparse -----------
disp('  ')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ')
disp(' function AlgParm=LaheBeamDFI(baasi0,Ntoerkts,esQkoormus,esFjoud,sSolmF,tsolm,tSiire,krdn,selem) ')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ')
disp('  ')

AlgParm=LaheBeamDFI(baasi0,Ntoerkts,esQkoormus,esFjoud,sSolmF,tsolm,tSiire,krdn,selem);
disp('  ')
AlgPar=AlgParm;
disp('  ')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ')
disp(' The element displacements and forces  determinated by the transfer matrix ')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ')
disp('  ')
#-------------The element displacements and forces----------
#
  suurused=[' displacement w - ';
            ' rotation    fi - ';
            ' shear force  Q - ';
            ' moment force M - '];
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for i=1:NEARV
 qkoormus=zeros(1,3);
Fjoud=zeros(1,2);
spvF=zeros(4,8);
vB=zeros(4,1);
vFz=zeros(4,1);
#
krda=i;
vF=zeros(4,8);
    EI=selem(i,9);   #  from topology
    GAr=selem(i,10); #    "      "  
Li=lvarras(i,1);
#
#AL#qz=qzZ(i,1);
#AL#aLx=aLXx(i,1);
#AL#Fz=FZz(i,1);
#
Fjoud=esFjoud(:,1:2,i);
qkoormus=esQkoormus(:,1:3,i);
xsamm=Li/Nmitmeks;  # element is divided into 4-th parts
xx=0;
AP=AlgPar(i,:)';
#---------The transfer matrix equation --------
for ij=1:Nmitmeks+1  # 5 - displacements and forces at x=0.0
##vvF=ylfmhvII(1.0,SII,xx,Li,EA,GAr,EI);
vvF=ylTfhlin(1.0,xx,GAr,EI);               # ylfhlin(1.0,xx,EA,GAr,EI);
#A#vvB=yzThqz(1.0,xx,qz,EI);               # yzhqz(1.0,xx,qx,qz,EA,EI);
##disp(' koondatud jõud ')
#A#vvFz=yzTfzv(1.0,xx,aLx,Fz,EI);          # yzfzv(1.0,xx,aLx,Fx,Fz,EA,EI);
##disp('jaotatud + koondatud jõud ')
vvB=ESTtalaKrmus(1.0,xx,Li,Fjoud,qkoormus,EI);
#A#vvB=vvB+vvFz;
Fvv(:,ij)=vvF*AP+vvB;
xx=xx+xsamm;
endfor
#%%%%%%%%%%%%%%%%%
# yzTfzv(baasi0,x,a,Fz,EJ)
# yzThqz(baasi0,x,qz,EJ))
# ylTfhlin(baasi0,x,GAr,EJ)
# yspTlvfmhvI(baasi0,x,l,GAr,EJ)
# yspTlfhlin(baasi0,x,GAr,EJ)
##%%%%%%%%%%%%%%%
#Fvv(1:2,:)
#Fvv(2:4,:)
VardaNr=i;
#disp(' Element Li is divided into 5-th parts ')
disp(sprintf('%15s %2i %17s %8.5f %5s', 'The displacements and forces of element nr',VardaNr,'  length of element ',Li,'m, '))
disp(sprintf('  element is divided into 5-th parts '))
#Li
#disp('  ')
#
for i=1:2 
disp(sprintf('%14s %9.3e  %9.3e  %9.3e  %9.3e  %9.3e  %9.3e',suurused(i,:), Fvv(i,1),   Fvv(i,2),   Fvv(i,3),  Fvv(i,4),  Fvv(i,5),  Fvv(i,6)))
endfor
#
for i=3:4 
disp(sprintf('%14s %9.3f  %9.3f  %9.3f  %9.3f  %9.3f  %9.3f',suurused(i,:), Fvv(i,1),   Fvv(i,2),   Fvv(i,3),  Fvv(i,4),  Fvv(i,5),  Fvv(i,6)))
endfor
 #
 #disp('------------------')
#
#for i=4:6 
#disp(sprintf('%16s %9.5f   %9.5f   %9.5f   %9.5f  %9.5f',suurusedNQ(i,:), FvvNQ(i,1),   %FvvNQ(i,2),   FvvNQ(i,3),  FvvNQ(i,4),  FvvNQ(i,5)))
#endfor
#
disp('----------------------------------------------------------------------------')
#%%%%%%%%%%%%%%
endfor
#
#%%%%%%%%
#SF=SisejoudPunktis(VardaNr,X,AlgPar,lvarras,selem,qxZ,qzZ,aLXx,FZz,FZx,suurused)  
SisejoudTalaPunktis(1,4.0-eps1,AlgPar,lvarras,selem,esFjoud,esQkoormus,suurused);
disp('-------------------------------------')
SisejoudTalaPunktis(2,4.0,AlgPar,lvarras,selem,esFjoud,esQkoormus,suurused);
disp('-------------------------------------')
SisejoudTalaPunktis(4,2.0,AlgPar,lvarras,selem,esFjoud,esQkoormus,suurused);
disp('-------------------------------------')
SisejoudTalaPunktis(6,4.0-eps1,AlgPar,lvarras,selem,esFjoud,esQkoormus,suurused);
disp('-------------------------------------')
SisejoudTalaPunktis(6,4.0,AlgPar,lvarras,selem,esFjoud,esQkoormus,suurused);
disp('-------------------------------------')
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
#
for i=1:NEARV
LkoordN=selem(i,11);
AkoordN=selem(i,12);
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);
VGRz(i,1)=0.0;
VGRz(i,2)=0.0;
endfor
#
##%%%%%%%%%%%
figure(1)
ax = gca ();
#axy=[ 0   0   1   0.55];
#joonepaksus=0.50000
#joonepaksus=0.75000
#joonepaksus=3
#set (ax, "outerposition", axy);
#set (ax, "outerposition", [ 0   0   1   0.55])
#ax = gca ();
#set (ax, "linewidth", joonepaksus);
#
grid off
#
Lpikk=Lp;
Hkrg=10;
#d=4.0;
d=4.0;
#
for j=1:NEARV
       IR=j;
 # 
#   axis([-d Lpikk+d -H-2.0 2.0], "ij")
   axis([-d Lpikk+d  -Hkrg d], "ij")
 #
plot(VGRx(IR,:),VGRz(IR,:),"3")
hold on
endfor
#
jaotT=2;
#
title('spESTGerberBeamWFi',"fontsize",12)
ax=gca();
      xlabel('x',"fontsize",12)
      ylabel('z',"fontsize",12)
 xticks(1,:)=-d:jaotT:Lpikk+d;
 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), 0.0+0.25, TST)
%text(krdn(IR,1), krdn(IR,2)+0.25, 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, 0.0, TST)
#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, 'Talaskeem')
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xvalg=5.0;
yvalg=-7.5;
sammuga=0.7;
text(5.0, yvalg-1.5*sammuga, 'The displacements and the forces numbers (daf''s) ')
#TSA=[TSA1 Tyhk TSA2 Tyhk TSA3] ; 
text(2.5, yvalg, 'w fi Q M at the beginning')
##xvalg=xvalg+sammuga;
yvalg=yvalg+sammuga;
for i=1:NEARV;     
#siireVardaA  joudVardaA(i,1:3)
TTA1=eval(sprintf('%3i',siireVardaA(i,1)));
TTA2=eval(sprintf('%3i',siireVardaA(i,2)));
TTA3=eval(sprintf('%3i',joudVardaA(i,1)));
TTA4=eval(sprintf('%3i',joudVardaA(i,2)));
   TSA1=int2str(TTA1);
   TSA2=int2str(TTA2);   
   TSA3=int2str(TTA3);  
   TSA4=int2str(TTA4);   
   Tyhk='  ';
TSA=[TSA1 Tyhk TSA2 Tyhk TSA3 Tyhk TSA4] ;                                                          
text(xvalg, yvalg, TSA)
##xvalg=xvalg+sammuga;
yvalg=yvalg+sammuga;
endfor
yvalg=yvalg+sammuga;
text(5.0, yvalg, 'The support reactions: 57 58 59 60 61  ')
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xvalg=17.0;
yvalg=-7.5;
sammuga=0.7;
#TSA=[TSA1 Tyhk TSA2 Tyhk TSA3] ; 
text(15.5, yvalg, 'w fi Q M  at the end')
##xvalg=xvalg+sammuga;
yvalg=yvalg+sammuga;
for i=1:NEARV;
#siireVardaL  joudVardaL(i,1:3)
TTA1=eval(sprintf('%3i',siireVardaL(i,1)));
TTA2=eval(sprintf('%3i',siireVardaL(i,2)));
TTA3=eval(sprintf('%3i',joudVardaL(i,1)));
TTA4=eval(sprintf('%3i',joudVardaL(i,2)));
   TSA1=int2str(TTA1);
   TSA2=int2str(TTA2);   
   TSA3=int2str(TTA3);  
   TSA4=int2str(TTA4);  
   Tyhk='  ';
TSA=[TSA1 Tyhk TSA2 Tyhk TSA3 Tyhk TSA4] ; 
text(xvalg, yvalg, TSA)
##xvalg=xvalg+sammuga;
yvalg=yvalg+sammuga;
endfor
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#
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','spESTGerberBeamWFi.eps');
#print('-dpng','spESTGerberBeamWFi.png');
print('-dfig','-landscape','-mono','-solid','spESTGerberBeamWFi.fig');
refresh
#
##%%%%%%%%%%%%%%%%%%%%%%%%%%
#spA
figure(2)
spy(spA)
title('spESTGerberTalaWFi spy(spA) - the sparse matrix spA(61,61)  non zero elements [4.1%] ')
refresh
#
#print('-deps','spESTtala_h6remaatriks.eps');
#print('-dpng','spESTtala_h6remaatriks.png');
print('-dfig','-landscape','-mono','-solid','spESTGerberTalaWFiaLahe_h6remaatriks.fig');
#refresh
#
# octave
#disp('-----------------------------------')
disp('  ')
disp(' The END ')
#%%%%%%%%%%%%%%%
 
