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.
45 lines
1.3 KiB
45 lines
1.3 KiB
%--------------------------
|
|
% @Author: Jingqiao Hu
|
|
% @Date: 2022-04-14 17:43:46
|
|
% @LastEditTime: 2022-04-25 10:42:27
|
|
%--------------------------
|
|
function Nh = init_Nh_triangle(nodes, faces)
|
|
x = nodes(:, 1); % [n, 1]
|
|
x = x(faces); % [m, 3]
|
|
y = nodes(:, 2); % [n, 1]
|
|
y = y(faces); % [m, 3]
|
|
[a, b, c, area] = paras(x, y); % [n,3]
|
|
|
|
Nh = zeros(size(x,1), 3, 3); % [m,3,3], [ele_num, nodes, 3 for 6]
|
|
for i = 1:3
|
|
% Nhj = zeros(size(x,1),1,3);
|
|
for j = 1:3
|
|
Nh(:,i,j) = (a(:,j) + b(:,j) .* x(:,i) + c(:,j) .* y(:,i)) / 2 ./ area;
|
|
end
|
|
% Nh(:,i,:) = kron(Nhj, eye(2)); % [2,6]
|
|
end
|
|
% Nh1 = kron(Nh, eye(2)); % [m,6]
|
|
end
|
|
|
|
function [a,b,c, area] = paras(x, y)
|
|
% a = [a1, a2, a3], n*3 , ...
|
|
a1 = x(:,2) .* y(:,3) - x(:,3) .* y(:,2);
|
|
a2 = x(:,3) .* y(:,1) - x(:,1) .* y(:,3);
|
|
a3 = x(:,1) .* y(:,2) - x(:,2) .* y(:,1);
|
|
a = [a1, a2, a3];
|
|
|
|
b1 = y(:,2) - y(:,3);
|
|
b2 = y(:,3) - y(:,1);
|
|
b3 = y(:,1) - y(:,2);
|
|
b = [b1, b2, b3];
|
|
|
|
c1 = -x(:,2) + x(:,3);
|
|
c2 = -x(:,3) + x(:,1);
|
|
c3 = -x(:,1) + x(:,2);
|
|
c = [c1, c2, c3];
|
|
|
|
area = 0.5 * ((x(:,1).* y(:,2) - x(:,2).*y(:,1)) + ...
|
|
(x(:,2).* y(:,3) - x(:,3).*y(:,2)) + ...
|
|
(x(:,3).* y(:,1) - x(:,1).*y(:,3)));
|
|
area = abs(area); % important!
|
|
end
|