##Naide4_6Uforce.m
#  vt lk 353
close all 
clear
## 
#dbwhere(750)
#
##=========================================================================
## This Program is written by Andres Lahe,   2015-03-30
##                    e-mail: andres.lahe@ttu.ee
## LAST MODIFIED: Andres Lahe,    2015-07-13
## Copyright (c)  2015  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
##=========================================================================
#
## lk 383 vt  ml/2=ql/2   m=q
## lk 343 bimomendi epüüri peallekandmint
##lk392 deplanatsiooni märklk398
 eps1=1.0*10e-11;
 eps2=0.00001 # 4 nulli pärast koma, eristab sortimisel  299.99999
 global spA
 global X
 global AP 
## global Fvvs
NNK=48 %%
#--------- 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('  Bvb=zeros(NNK,1);  ')
Bvb=zeros(NNK,1);   # nullistan võrrandisüsteemi vabaliikmed
disp('   ')
#######
#baasi0=1.0   # baasjäikus
baasi0=1.0e+8  # baasjäikus
l1=400  # cm   # varda 1 pikkus   
l2=500  # cm   # varda 2 pikkus 
l2p=l2/2  # cm   # varda 2 pool pikkust 
l3=400  # cm   # varda 1 pikkus 
#x=0.0 # cm   # koordinaat xs
#
qz=10.0
mx=10.0 # N    # lausmoment 
GAr=  2.8348e+13
#
GIt1= 1.5637e+09
EIw1= 2.8348e+13 #Ncm4 # kooldumisjäikus
k1=sqrt(GIt1/EIw1) ## k =  0.0074270
EI1 = EIw1/16.097

GIt2= 1.7753e+09
EIw2= 2.9254e+13 #Ncm4 # kooldumisjäikus
k2=sqrt(GIt2/EIw2) ## k =  0.007790
EI2 = EIw2/16.097
######################################
#i1w=EIw1/l1
#i1x=EI1/l1
#i2w=EIw2/l2p
#i2x=EI2/l2p
#####################################
#qz=100.0
#mx=100.0 # N    # lausmoment 
#mx1=mx;EIX =  1.7611e+12
#mx=1000.0 
#mx=10000.0 # N    # lausmoment
#Mx2=3.20e+05 # Ncm    # koondmoment 320 000
#Mx2=3.20e+05 # Ncm  -147 m kgcm  lk 352
#Mx=3.20e+05 # Ncm    # koondmoment 
##By=-1.00e+07 # Ncm2    # koondbimoment 10 000 000
#By=-1.00e+07
##By=1.00e+11
#a2=300  # cm   # lausmomendi algus
# Iw=1349900
# It= 195.50
#GIt= 1.5837e+09  # Ncm2      # väändejäikus
#######
#G =  8000000                                                                                                                                                                                                            
#E =  21000000                                                                                                                                                                                                             
#It =  195.50                                                                                                 
#Iw =  1349900                                                                                                       
#GIt=G*It =  1.5640e+09 
#EIw=E*Iw =  2.8348e+13                                                                                                
#kt=sqrt(GIt/EIw)   =  0.0074270
#k=sqrt(GIt/EIw)
#######
#Iw=1349900  # cm6
#Ix=83860  # cm4
##EIw=2.8356e+13                                                                                                                                                                                                                                                                                                                  
#EIX=E*Ix
#EIX =  1.7611e+12
#EI =  1.7611e+12
#EJw=E*Iw
#EJw =  2.8348e+13
##################################################################
#==========
#    Nodal coordinates
#==========
krdn=[#  x       y
        0.0     0.0;    # node 1
        l1      0.0;    # node 2
        l1      -l2;    # node 3        
        0.0     -l2];   # node 4 
#==========
#%%%%%%%%%%%%
SARV=size(krdn);
NSARV=SARV(1,1);
#%%%
for i=1:NSARV
sJrN(i)=i;
end
#
sJrNT=sJrN';
###%%#i1w=EIw1/l1
#i1x=EI1/l1
#i2w=EIw2/l2p
#i2x=EI2/l2p
#==========
selemjl=[# n1 - beginning of the element,
#              n2 - end of the element
           1   2;    # element 1
           2   3;    # element 2 
           3   4];   # element 3  
#==========
#%%%%%%%%%%%%
selem=[selemjl(:,1) selemjl(:,2)];
EARV=size(selem);
NEARV=EARV(1,1);
#%%%%%%%%%%%%
selem=[selemjl(:,1) selemjl(:,2)];
##%%%
for i=1:NEARV
JrN(i)=i;
end
%
JrNT=JrN';
#%%%%
disp('===============================')
disp(' Nodal coordinates ')
disp('JrNr        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('-----------------------------')
#%%%
##%%%
disp('=========================')
disp(' Topology of elements    ')
disp('     At the     At the   ')
disp(' Nr  beginning, end  ')
disp('-------------------------')
#disp([selem])
for i=1:NEARV
disp(sprintf('  %2i     %2i    %2i',JrNT(i),    selemjl(i,1), selemjl(i,2)))
endfor
disp('-------------------------')
#
###
for i=1:NEARV
LkoordN=selem(i,2);
AkoordN=selem(i,1);
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
#
#Varda pikkus
#lvarras
for i=1:NEARV
lvarras(i,1)=sqrt(DeltaX(i)^2+DeltaZ(i)^2);
endfor
#
for i=1:NEARV
cosAlpha(i,1)=DeltaX(i)/lvarras(i,1);
cosBeta(i,1)=DeltaZ(i)/lvarras(i,1);
%
suunakosin= zeros(2,2);
        suunakosin(1,1)=cosAlpha(i,1);
        suunakosin(1,2)=-cosBeta(i,1);
                
        suunakosin(2,1)=cosBeta(i,1);  
        suunakosin(2,2)=cosAlpha(i,1);  
                    
%suunakosin
%TJ=suunakosin(NEARV,cosAlpha(i,1),cosBeta(i,2));
TM2x2(:,:,i)=suunakosin;
endfor
%
disp(' Varda 1 teisendus maatriks ')
TMI2x2=TM2x2(:,:,1)
disp(' Varda 2 teisendus maatriks ')
TMII2x2=TM2x2(:,:,2)
disp(' Varda 3 teisendus maatriks ')
TMIII2x2=TM2x2(:,:,3)
#
#disp('-----------------------------')
disp('===================================')
disp(' The direction cosines of elements ')
disp('Element  cosAlpha  cosBeta ')
disp('-----------------------------------')
#disp([cosAlpha cosBeta])
for i=1:NEARV
disp(sprintf('  %2i     %7.4f   %7.4f',JrNT(i),    cosAlpha(i,:), cosBeta(i,:)))
endfor
disp('-----------------------------------')
#
#
%
%%%%%%%%
#
i1w=EIw1/l1
#i1w=1.0
i1x=EI1/l1
 disp(' ')
i2w=EIw2/l2p
i2x=EI2/l2p
 disp(' ')
k1=sqrt(GIt1/EIw1)
k2=sqrt(GIt2/EIw2)
% 
nu1=l1*sqrt(GIt1/EIw1)
 ch1=cosh(nu1);
 sh1=sinh(nu1);
 th1=tanh(nu1)
%%%%%%%%
 nu2=l2p*sqrt(GIt2/EIw2)
%
 ch2=cosh(nu2);
 sh2=sinh(nu2);
 th2=tanh(nu2) 
%%%%%%%%%%%%%%%%% lk457 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 disp(' ')
lamb1=(sh1*nu1^3)/(nu1*sh1-2*ch1+2)
lamb2=(sh2*nu2^3)/(nu2*sh2-2*ch2+2)
b1=th1/nu1
b2=th2/nu2
mu1=(nu1*(nu1*ch1-sh1))/(nu1*sh1-2*ch1+2)
        mu1=i1w*mu1
mu2=(nu2*(nu2*ch2-sh2))/(nu2*sh2-2*ch2+2)
alp1=((ch1-1)*nu1^2)/(nu1*sh1-2*ch1+2)
alp1=i1w*alp1
alp2=((ch2-1)*nu2^2)/(nu2*sh2-2*ch2+2)
r1=(nu1*ch1-sh1)/(nu1^2*sh1)
r2=(nu2*ch2-sh2)/(nu2^2*sh2)
s1=(sh1-nu1)/(nu1^2*sh1)
s2=(sh2-nu2)/(nu2^2*sh2)
f1=(ch1-1)/(nu1*sh1)
f2=(ch2-1)/(nu2*sh2)
mu1=(nu1*(nu1*ch1-sh1))/(nu1*sh1-2*ch1+2)
        mu1=i1w*mu1
a1=1/ch1
a1=1/ch2
gamm1=(nu1*(sh1-nu1))/(nu1*sh1-2*ch1+2)
        gamm1=i1w*gamm1
gamm2=(nu2*(sh2-nu2))/(nu2*sh2-2*ch2+2)

 disp(' ')
%%%%%%%%% 
D= lamb1*(b2+i2w*b1)+l1^2*i2x*b1*(mu1*b2+i2w) 
H=D/2 
%
 disp(' ')
Mxc=-mx*l2p^2*l1*i2x*(alp1*b1*b2*r2)/D 
Mxd=-mx*l2p
 
  disp(' ')
Bc = -mx*l2p^2*b2*(lamb1*(s2+i2w*f2*b1)+ l1^2*i2x*b1*(mu1*s2+i2w*f2))/D  
Bda = mx*l2p^2*b2*r2*(lamb1+l1^2*i2x*mu1*b1)/D
Beb = Bda
Bad = mx*l2p^2*b2*r2*(a1*lamb1-l1^2*i2x*b1*gamm1)/D   
Bbe = Bad 
 
%%%%%%%%%  
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 