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.

91 lines
2.8 KiB

3 years ago
%--------------------------
% @Author: Jingqiao Hu
% @Date: 2021-11-26 21:02:20
% @LastEditTime: 2022-01-03 23:19:44
% get the bnid node idx for each cvt
% NOTE: the bdofs must CONFORMS to phi
%--------------------------
function [bdofs_cell, dofid_cell, nbnodes_cell] = border_check(c_faces, c_nodes, c_vedges)
cvt_num = size(c_faces, 1);
bdofs_cell = cell(cvt_num, 1); % saves each bdofs
dofid_cell = cell(cvt_num, 1);
nbnodes_cell = cell(cvt_num, 1); % saves the number of boundary dofs
% figure;
% hold on;
for ele = 1 : cvt_num
nodes = c_nodes{ele};
vedges = c_vedges{ele};
bnid = [];
edge_num = size(vedges, 1);
bnum = zeros(edge_num, 1);
% NOTE: DO NOT using parfor! sequence related
for i = 1 : edge_num
n1 = vedges(i, 1:2);
n2 = vedges(i, 3:4);
bnid_each = [];
for j = 1 : size(nodes, 1)
p = nodes(j, :);
flag = point_on_seg(n1, n2, p);
if flag == 1
bnum(i) = bnum(i) + 1;
bnid_each = [bnid_each; j];
end
end
bnid_new = rearrange_border_nodes(n1, nodes, bnid_each);
bnid = [bnid; bnid_new(1:end-1)]; % remove the adjoint corner node
% bnid = [bnid; bnid_each(1:end-1)]; % remove the adjoint corner node
end
nbnodes_cell{ele} = bnum;
% figure(1); clf;
% scatter(nodes(bnid, 1), nodes(bnid, 2), 'filled'); hold on;
% faces = c_faces{ele};
% trimesh(faces, nodes(:,1), nodes(:,2));
bdofs = [2*bnid-1, 2*bnid]';
bdofs_cell{ele} = bdofs(:);
alldofs = 1 : size(nodes, 1) * 2;
idofs = setdiff(alldofs, bdofs(:));
newdofs = [bdofs(:)', idofs];
dofid = reorder_dofs(newdofs);
dofid_cell{ele} = dofid;
end
end
% rearrange as the distance to node n1
function bnid_new = rearrange_border_nodes(n1, nodes, bnid)
nodes0 = nodes(bnid, :);
dist = (nodes0(:,1) - n1(1)).^2 + (nodes0(:,2) - n1(2)).^2;
nodes0(:,3) = dist;
nodes0(:,4) = bnid;
nodes1 = sortrows(nodes0, 3);
bnid_new = nodes1(:, 4);
end
% segment = [n1, n2]
% if point p on the segment
function flag = point_on_seg(n1, n2, p)
delta = 1e-2;
if abs((p(1) - n1(1)) * (n2(2) - n1(2)) - (n2(1) - n1(1)) * (p(2) - n1(2))) < delta... % cross product
&& (min(n1(1), n2(1)) < p(1) || abs(min(n1(1), n2(1)) - p(1)) < delta)...
&& (p(1) < max(n1(1), n2(1)) || abs(p(1) - max(n1(1), n2(1))) < delta) ...
&& (min(n1(2), n2(2)) < p(2) || abs(min(n1(2), n2(2)) - p(2)) < delta) ...
&& (p(2) < max(n1(2), n2(2)) || abs(p(2) - max(n1(2), n2(2))) < delta)
flag = 1;
else
flag = 0;
end
end