a 2D version
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 

101 lines
3.0 KiB

%--------------------------
% @Author: Jingqiao Hu
% @Date: 2020-10-29 15:04:15
% @LastEditTime: 2021-04-19 10:11:20
%--------------------------
% gradient of shape function on one ele, four nodes, N1 ~ N4
%
% gradN: [4*8]
% [dN1dx, 0, dN2dx, 0, ...
% 0, dN1dx, 0, dN2dx, ...
% dN1dy, 0, dN2dy, 0, ...
% 0, dN1dy, 0, dN2dy, ...]
% This sequence is setted for computaion
%
% gradN2: [3*8]
% [dN1dx, 0, dN2dx, 0, ...
% 0, dN1dy, 0, dN2dy, ...
% dN1dy, dN1dx, dN2dy, dN2dx]
function [gradN, gradN2, N] = shape_gradient_ele(a, b, sample_num)
if sample_num==1
[gradN] = shape_gradient_1(a, b);
N = [];
gradN2 = [];
else
[gradN, gradN2, N] = shape_gradient_4(a, b);
end
end
% 4*8
function [gradN] = shape_gradient_1(a, b)
GN_x = 0;
GN_y = 0;
dN_x = 1/4*[-(1-GN_x) (1-GN_x) (1+GN_x) -(1+GN_x)];
dN_y = 1/4*[-(1-GN_y) -(1+GN_y) (1+GN_y) (1-GN_y)];
J = [dN_x; dN_y]*[ -a a a -a;
-b -b b b]';
G = [inv(J) zeros(size(J));
zeros(size(J)) inv(J)];
dN(1,1:2:8) = dN_x; dN(3,1:2:8) = dN_y;
dN(2,2:2:8) = dN_x; dN(4,2:2:8) = dN_y;
B = G*dN;
gradN = B;
end
% 4 * 8 & 3*8 for every cell
function [gradN, gradN2, N] = shape_gradient_4(a, b)
GaussNodes = [-1/sqrt(3); 1/sqrt(3)];
gradN = cell(4,1);
gradN2 = cell(4,1);
L = [1 0 0 0; 0 0 0 1; 0 1 1 0];
nodes = [-1, 1];
N = cell(4,1);
for i = 1:2
for j = 1:2
gp = 2*(i-1)+j;
GN_x = GaussNodes(i);
GN_y = GaussNodes(j);
% GN_x = nodes(i); GN_y = nodes(j);
dN_x = 1/4*[-(1-GN_x) (1-GN_x) (1+GN_x) -(1+GN_x)];
dN_y = 1/4*[-(1-GN_y) -(1+GN_y) (1+GN_y) (1-GN_y)];
J = [dN_x; dN_y]*[ -a a a -a;
-b -b b b]';
G = [inv(J) zeros(size(J));
zeros(size(J)) inv(J)];
dN1(1,1:2:8) = dN_x; dN1(2,1:2:8) = dN_y;
dN1(3,2:2:8) = dN_x; dN1(4,2:2:8) = dN_y;
% dNdx = 4*8
B = G*dN1;
gradN{gp} = B([1,3,2,4],:);
% dNdx = 3*8
dN1(1,1:2:8) = dN_x; dN1(2,1:2:8) = dN_y;
dN1(3,2:2:8) = dN_x; dN1(4,2:2:8) = dN_y;
B2 = L*G*dN1;
gradN2{gp} = B2;
% N = 2*8
GN_x = nodes(i); GN_y = nodes(j);
N1 = (1-GN_x)*(1-GN_y)/4;
N2 = (1+GN_x)*(1-GN_y)/4;
N3 = (1+GN_x)*(1+GN_y)/4;
N4 = (1-GN_x)*(1+GN_y)/4;
N{gp} = [N1, 0, N2, 0, N3, 0, N4, 0;
0, N1, 0, N2, 0, N3, 0, N4];
% N2{gp} = [N1, N1;
% N2, N2;
% N3, N3;
% N4, N4;]; % 4*2
end
end
end