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.

90 lines
2.5 KiB

3 years ago
%--------------------------
% @Author: Jingqiao Hu
% @Date: 2021-11-26 18:33:42
% @LastEditTime: 2022-04-29 15:27:29
%--------------------------
function [psi_cell] = preparation_cvt(e1, e0, penal, nbnodes_cell, phi_cell, dofid_cell, ...
rho_cell, faces_cell, c_nodes, sK_cell, c_edofMat, cbns)
ncvt = size(rho_cell, 1);
psi_cell = cell(ncvt, 1);
c_NH = cell(ncvt, 1);
% plot ele = 44
parfor ele = 1 : ncvt
faces = faces_cell{ele};
nele = size(faces, 1);
[iK, jK, edofMat, ~, ~] = forAssemble_trimesh(faces);
order1 = reshape(edofMat', [], 1);
rho = rho_cell{ele}(:); % rho(:) = 1;
phi = phi_cell{ele};
sK0 = sK_cell{ele};
dofid = dofid_cell{ele};
nedges = length(nbnodes_cell{ele});
nbdofs = (sum(nbnodes_cell{ele}) - nedges) * 2;
k = assemble_micro_k_hetero(dofid, rho, e0, e1, penal, sK0, iK, jK);
Me = build_M_simply(nbdofs, k, phi);
Re = [phi; Me]; % [2n, 2v]
R1 = Re(dofid(order1), :); % [6m, 2v]
R2 = reshape(R1, 6, nele, []); % [6, m, 2v]
R3 = permute(R2, [1,3,2]);
% edof = c_edofMat{ele};
% u = reshape(cbns', [], 1);
% ue = u(edof);
% umi = Re * ue;
% fineNodes = reshape(umi(dofid), 2, [])';
% figure; scatter(fineNodes(:,1), fineNodes(:,2)); hold on;
% re = reshape(ue, 2, [])';
% scatter(re(:,1), re(:,2), 'filled');
%
% dt = triangulation(faces, fineNodes);
% triplot(dt);
psi_cell{ele} = R3;
%% surface of NH
% nodes = c_nodes{ele};
% Nh = init_Nh_triangle(nodes, faces); % [m,3]
%
% n_vdofs = size(Re, 2);
%
% R3 = Re(dofid, :); % [2n, 2v]
% NHe = NH_alpha(edofMat, n_vdofs/2, Nh, R3, size(nodes, 1));
% NH{ele, 1} = NHe; % [n,2,2,v]
%
% plot_NH_cvt(nodes, faces, NHe);
end
end
function plot_NH_cvt(nodes, faces, NH)
% [n,2,2,v]
x = nodes(:, 1); % [n, 1]
% x = x(faces); % [m, 3]
y = nodes(:, 2); % [n, 1]
% y = y(faces); % [m, 3]
for i = 1 :size(NH,4)
figure;
% for j = 1:3
cv = NH(:,2,2,i) * 10;
surfir(x,y,cv,0);
% scatter3(x, y, cv);
% colormap(jet); caxis([0,1]);
% s.FaceColor = 'w';
% s.EdgeColor = [0.2,0.2,0.2];
% hold on;
grid off; axis equal; axis off;
set(gca,'FontSize',14,'FontWeight','bold');
% end
end
end