% 接続行列の計算例

a = [ -1;  1;  0 ];
b = [ -1;  0;  1 ];
h = 2;
S = 1/2;

luu = (h*S)*a*a'; luv = (h*S)*a*b';
lvu = (h*S)*b*a'; lvv = (h*S)*b*b';
muu = 2*luu + lvv; muv = lvu;
mvu = luv; mvv = 2*lvv + luu;

fprintf("行列 H_lambda^ABC, H_mu^ABC\n");
l = [ luu, luv; lvu, lvv ];
m = [ muu, muv; mvu, mvv ];
l
m

fprintf("部分行列 J_lambda^ABC, J_mu^ABC\n");
Jlambda_ijk = l([1,4,2,5,3,6], [1,4,2,5,3,6]);
Jmu_ijk     = m([1,4,2,5,3,6], [1,4,2,5,3,6]);
Jlambda_ijk
Jmu_ijk

triangle = [ 1, 2, 4; ...
             5, 4, 2; ...
             2, 3, 5; ...
             6, 5, 3 ];

Jlambda = zeros(2*6, 2*6);
Jmu = zeros(2*6, 2*6);

for row = triangle'
    i = row(1); j = row(2); k = row(3);
    %Jlambda(2*i-1:2*i,2*i-1:2*i) = Jlambda(2*i-1:2*i,2*i-1:2*i) + Jlambda_ijk(1:2,1:2);
    %Jlambda(2*i-1:2*i,2*j-1:2*j) = Jlambda(2*i-1:2*i,2*j-1:2*j) + Jlambda_ijk(1:2,3:4);
    %Jlambda(2*i-1:2*i,2*k-1:2*k) = Jlambda(2*i-1:2*i,2*k-1:2*k) + Jlambda_ijk(1:2,5:6);
    %Jlambda(2*j-1:2*j,2*i-1:2*i) = Jlambda(2*j-1:2*j,2*i-1:2*i) + Jlambda_ijk(3:4,1:2);
    %Jlambda(2*j-1:2*j,2*j-1:2*j) = Jlambda(2*j-1:2*j,2*j-1:2*j) + Jlambda_ijk(3:4,3:4);
    %Jlambda(2*j-1:2*j,2*k-1:2*k) = Jlambda(2*j-1:2*j,2*k-1:2*k) + Jlambda_ijk(3:4,5:6);
    %Jlambda(2*k-1:2*k,2*i-1:2*i) = Jlambda(2*k-1:2*k,2*i-1:2*i) + Jlambda_ijk(5:6,1:2);
    %Jlambda(2*k-1:2*k,2*j-1:2*j) = Jlambda(2*k-1:2*k,2*j-1:2*j) + Jlambda_ijk(5:6,3:4);
    %Jlambda(2*k-1:2*k,2*k-1:2*k) = Jlambda(2*k-1:2*k,2*k-1:2*k) + Jlambda_ijk(5:6,5:6); 
    % 上記の9行をまとめて書く
    p = [2*i-1:2*i, 2*j-1:2*j, 2*k-1:2*k];
    Jlambda(p,p) = Jlambda(p,p) + Jlambda_ijk;
    Jmu(p,p) = Jmu(p,p) + Jmu_ijk;
end

fprintf("行列 J_lambda, J_mu\n");
Jlambda
Jmu
