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.
95 lines
2.9 KiB
95 lines
2.9 KiB
%--------------------------
|
|
% @Author: Jingqiao Hu
|
|
% @Date: 2021-11-26 21:02:20
|
|
% @LastEditTime: 2022-04-28 10:47:50
|
|
|
|
% for each cvt and each surface: get the boundary node idx
|
|
% NOTE: the bdofs must CONFORMS to phi
|
|
% CHECK: bdofs order, do we need rearrange?
|
|
% paras:
|
|
% c_polygon: all segs for each face and each cell
|
|
%--------------------------
|
|
function [c_bnid, c_dofid, c_bnidSaved, c_bnidSavedID, nbnidFine] = border_check(c_node, c_polygon)
|
|
|
|
ncell = size(c_node, 1);
|
|
c_dofid = cell(ncell, 1);
|
|
c_bnid = cell(ncell, 1);
|
|
c_bnidSavedID = cell(ncell, 1);
|
|
c_bnidSaved = cell(ncell, 1);
|
|
nbnidFine = zeros(ncell, 1);
|
|
|
|
% figure;
|
|
for ele = 1 : ncell
|
|
nodes = c_node{ele};
|
|
polys = c_polygon{ele};
|
|
|
|
[c_bnid_face, c_bnidSaved_face, c_bnidSavedID_face, bnid] = boundary_nodes(nodes, polys);
|
|
|
|
c_bnid{ele} = c_bnid_face;
|
|
c_bnidSavedID{ele} = c_bnidSavedID_face;
|
|
c_bnidSaved{ele} = c_bnidSaved_face;
|
|
|
|
nbnidFine(ele) = length(bnid);
|
|
% bnid = sort(bnid); % important, phi-B is order as this
|
|
bdofs = [3*bnid-2, 3*bnid-1, 3*bnid]';
|
|
|
|
alldofs = 1 : size(nodes, 1) * 3;
|
|
idofs = setdiff(alldofs, bdofs(:));
|
|
newdofs = [bdofs(:)', idofs];
|
|
dofid = reorder_dofs(newdofs);
|
|
c_dofid{ele} = dofid;
|
|
end
|
|
end
|
|
|
|
function [c_bnid, c_bnidSaved, c_bnidSavedID, bnid] = boundary_nodes(nodes, c_poly)
|
|
|
|
tol = 1e-3;
|
|
|
|
nface = size(c_poly, 1);
|
|
npnt = size(nodes, 1);
|
|
|
|
% bnid for each face
|
|
c_bnid = cell(nface, 1);
|
|
% bnid for each face, removed duplicated node with previous
|
|
c_bnidSaved = cell(nface, 1);
|
|
% index of bnidSaved, for phi-B
|
|
c_bnidSavedID = cell(nface, 1);
|
|
bnid = [];
|
|
|
|
% figure; scatter3(nodes(:,1), nodes(:,2), nodes(:,3),'b'); hold on;
|
|
% NOTE: DO NOT using parfor! sequence related
|
|
for i = 1 : nface
|
|
segs = c_poly{i};
|
|
|
|
p1 = segs(1, 1:3);
|
|
v1 = segs(1, [4:6]) - segs(1, [1:3]);
|
|
v2 = segs(2, [4:6]) - segs(2, [1:3]);
|
|
|
|
% TODO: change to surface
|
|
|
|
plane = repmat([p1, v1, v2], npnt, 1);
|
|
projNode = projPointOnPlane(nodes, plane);
|
|
|
|
dist = projNode - nodes;
|
|
[idx, ~] = find(abs(dist(:, 1)) < tol & ...
|
|
abs(dist(:, 2)) < tol & ...
|
|
abs(dist(:, 3)) < tol);
|
|
|
|
c_bnid{i} = idx; % bnid of this face
|
|
|
|
[bnidSaved, ia] = setdiff(idx, bnid, 'stable');
|
|
bnid = [bnid; bnidSaved(:)];
|
|
|
|
% id of bnidSaved of this face
|
|
% [~,bnidSavedIDX,~] = intersect(idx, bnidSaved, 'rows', 'stable');
|
|
|
|
c_bnidSavedID{i} = ia;
|
|
c_bnidSaved{i} = length(bnid) - length(bnidSaved) + 1 : length(bnid);
|
|
|
|
% x = segs(:,[1,4]);
|
|
% y = segs(:,[2,5]);
|
|
% z = segs(:,[3,6]);
|
|
% plot3(x',y',z','LineWidth',2); hold on;
|
|
end
|
|
% scatter3(nodes(bnid,1), nodes(bnid,2), nodes(bnid,3), 'filled', 'r'); hold on;
|
|
end
|