how to triple Integrate this?

1 view (last 30 days)
Saurabh Agarwal
Saurabh Agarwal on 7 Aug 2021
Answered: Walter Roberson on 8 Aug 2021
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;
for i = 1:24
for j = 1:24
K_mat(i,j) = @(xi,eta,zeta) K_ma(i,j);
K(i,j) = integral3(K_ma(i,j),-1,1,-1,1,-1,1);
end
end

Sign in to answer this question.

Accepted Answer

Walter Roberson
Walter Roberson on 8 Aug 2021
Ran in:
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)
K = 

More Answers (0)

Categories

Find more on Linear Least Squares in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!