% Mondai52.m
% 
% 5章の問題（2）更新枠法による要素力計算
% 
% Copyright (C) 2007 清水 信行, 曽我部 潔, 藤川 猛 All Rights Reserved.

% =================================================================================
% プログラムの使用条件
% 1. 本プログラムおよびデータを商用で，頒布してはならない．
% 2. 本プログラムおよびデータを著作者の許可なく営利目的で利用してはならない．
% 3. 本プログラムおよびデータの改変は，営利目的でない限り自由である．
% 4. 本プログラムおよびデータを使用することにより生じた損害などについては，著者，
%    プロクラム著作者，(社)日本機械学会，および(株)コロナ社は一切の責任を負わない．
% 5. 本プログラムの著作権は，本書編著者ならびにプログラム著作者に帰す．
% =================================================================================

com='Problem 5.2 Incremantal method'
L0=0.5*sqrt(2);EA=1.2e8;EI=4e3;                       % データ設定
Qi =[0.1; 0.1;pi/4]  ;Qj=[0.59;0.61;0.001+pi/4];      %現行座標
Qib=[0.1; 0.1;pi/4] ;Qjb=[0.60;0.60;pi/4];            %先行座標
Q=[Qi;Qj];Qb=[Qib;Qjb];
L =norm( Qj(1:2) -Qi(1:2));                           %現行長さ
Lb=norm(Qjb(1:2)-Qib(1:2));                            %先行長さ
Epsb=[0 0 0];                                          %先行変形量
co=cos(pi/4);si=sin(pi/4);Lamb=[co -si 0;si co 0; 0 0 1];
Tb=[Lamb zeros(3);zeros(3) Lamb];                      %先行変換マトリックス
    dQ=Q-Qb;                                           %全体枠変位増分
    dq=Tb*dQ;                                          %先行枠枠変位増分
    dw=(dq(5)-dq(2))/Lb
    du=(dq(4)-dq(1))/Lb
depsL=L-Lb;                                            %変形量増分
depsi=dq(3)-du*(1-dw);
depsj=dq(6)-du*(1-dw);
Eps=Epsb+[depsL depsi depsj]

fxi=6*EI/(L0*L)*(Eps(2)+Eps(3))                        %現行枠要素力
fzi=-EA/L0*Eps(1)
Mi= 2*EI/L0*(2*Eps(2)+Eps(3))
fxj=-6*EI/(L0*L)*(Eps(2)+Eps(3))
fzj=EA/L0*Eps(1)
Mj= 2*EI/L0*(Eps(2)+2*Eps(3))

dLam=[1-du^2/2 -du*(1-dw) 0;du*(1-dw) 1-du^2/2 0; 0 0 1];%
dT=[dLam zeros(3);zeros(3) dLam];                       ;%
T=dT*Tb;                                                ;%現行変換マトリックス
F=T'*[fxi;fzi;Mi;fxj;fzj;Mj]                             %現行要素力
