clear
clc
D = 200000;
syms xi eta zeta
N(1) = (1/8)*(1-xi)*(1-eta)*(1-zeta);
N(2) = (1/8)*(1+xi)*(1-eta)*(1-zeta);
N(3) = (1/8)*(1+xi)*(1+eta)*(1-zeta);
N(4) = (1/8)*(1-xi)*(1+eta)*(1-zeta);
N(5) = (1/8)*(1-xi)*(1-eta)*(1+zeta);
N(6) = (1/8)*(1+xi)*(1-eta)*(1+zeta);
N(7) = (1/8)*(1+xi)*(1+eta)*(1+zeta);
N(8) = (1/8)*(1-xi)*(1+eta)*(1+zeta);
for i = 1:8
N_xi(i) = diff(N(i),xi);
N_eta(i) = diff(N(i),eta);
N_zeta(i) = diff(N(i),zeta);
end
% Nodes xi eta zeta
% 1 -1 -1 -1
% 2 1 -1 -1
% 3 1 1 -1
% 4 -1 1 -1
% 5 -1 -1 1
% 6 1 -1 1
% 7 1 1 1
% 8 -1 1 1
% x1_xi = (N_xi(1)*-1) + (N_xi(2)*1) + (N_xi(3)*1) + (N_xi(4)*-1) + (N_xi(5)*-1) + (N_xi(6)*1) + (N_xi(7)*1) + (N_xi(8)*-1);
% x1_eta = (N_eta(1)*-1) + (N_eta(2)*1) + (N_eta(3)*1) + (N_eta(4)*-1) + (N_eta(5)*-1) + (N_eta(6)*1) + (N_eta(7)*1) + (N_eta(8)*-1);
% x1_zeta = (N_zeta(1)*-1) + (N_zeta(2)*1) + (N_zeta(3)*1) + (N_zeta(4)*-1) + (N_zeta(5)*-1) + (N_zeta(6)*1) + (N_zeta(7)*1) + (N_zeta(8)*-1);
% x2_xi = (N_xi(1)*-1) + (N_xi(2)*-1) + (N_xi(3)*1) + (N_xi(4)*1) + (N_xi(5)*-1) + (N_xi(6)*-1) + (N_xi(7)*1) + (N_xi(8)*1);
% x2_eta = (N_eta(1)*-1) + (N_eta(2)*-1) + (N_eta(3)*1) + (N_eta(4)*1) + (N_eta(5)*-1) + (N_eta(6)*-1) + (N_eta(7)*1) + (N_eta(8)*1);
% x2_zeta = (N_zeta(1)*-1) + (N_zeta(2)*-1) + (N_zeta(3)*1) + (N_zeta(4)*1) + (N_zeta(5)*-1) + (N_zeta(6)*-1) + (N_zeta(7)*1) + (N_zeta(8)*1);
% x3_xi = (N_xi(1)*-1) + (N_xi(2)*-1) + (N_xi(3)*-1) + (N_xi(4)*-1) + (N_xi(5)*1) + (N_xi(6)*1) + (N_xi(7)*1) + (N_xi(8)*1);
% x3_eta = (N_eta(1)*-1) + (N_eta(2)*-1) + (N_eta(3)*-1) + (N_eta(4)*-1) + (N_eta(5)*1) + (N_eta(6)*1) + (N_eta(7)*1) + (N_eta(8)*1);
% x3_zeta = (N_zeta(1)*-1) + (N_zeta(2)*-1) + (N_zeta(3)*-1) + (N_zeta(4)*-1) + (N_zeta(5)*1) + (N_zeta(6)*1) + (N_zeta(7)*1) + (N_zeta(8)*1);
%
% J = [x1_xi x1_eta x1_zeta;x2_xi x2_eta x2_zeta;x3_xi x3_eta x3_zeta];
J = [1 0 0;0 1 0;0 0 1];
Inv_J = inv(J);
Det_J = det(J);
N_x = [Inv_J(1,1) * N_xi; Inv_J(2,2) * N_eta; Inv_J(3,3) * N_zeta];
N_x = [N_x(1,1) 0 0 N_x(1,2) 0 0 N_x(1,3) 0 0 N_x(1,4) 0 0 N_x(1,5) 0 0 N_x(1,6) 0 0 N_x(1,7) 0 0 N_x(1,8) 0 0;...
0 N_x(2,1) 0 0 N_x(2,2) 0 0 N_x(2,3) 0 0 N_x(2,4) 0 0 N_x(2,5) 0 0 N_x(2,6) 0 0 N_x(2,7) 0 0 N_x(2,8) 0;...
0 0 N_x(3,1) 0 0 N_x(3,2) 0 0 N_x(3,3) 0 0 N_x(3,4) 0 0 N_x(3,5) 0 0 N_x(3,6) 0 0 N_x(3,7) 0 0 N_x(3,8)];
K_ma = N_x.' * N_x;
K = int(int(int(K_ma,xi,-1,1),eta,-1,1),zeta,-1,1)
