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.

89 lines
2.9 KiB

3 years ago
%--------------------------
% @Author: Jingqiao Hu
% @Date: 2021-11-27 15:21:49
% @LastEditTime: 2021-12-06 20:12:34
% could be arbitrary version without considering the length?
% NOTE: the bdofs must CONFORMS to phi
% bdofs order of phi: boundary dofs on each edge
% Three situations:
% 1. fine boundary nodes >= 2 on this edge, phi as usual
% 2. fine boundary node = 1 on this edge, phi with 1 contro node
% 3. fine boundary node = 0 on this edge, phi = eye(4)
%--------------------------
function phi_cell = parametric_B_cvt(nbnodes_cell, nbridge_each)
ncvt = size(nbnodes_cell, 1);
phi_cell = cell(ncvt,1);
% figure;
B2 = simple_quad_Bezier(3);
for ele = 1 : ncvt
nbele = nbnodes_cell{ele} - 1;
nedges = size(nbele, 1);
ncbn = 1;
nbdofs = 1;
standard_nedges = length(find(nbele>2));
ncbn_1v = standard_nedges * ((nbridge_each+1) * 2 + nbridge_each) + nedges;
ncbn_1v = ncbn_1v + sum(nbele(nbele <=2) - 1); % for situations 2 & 3
phi_1v = zeros(sum(nbele)*2, ncbn_1v * 2);
for edge = 1 : nedges
nb_edge = nbele(edge);
% three situations
if nb_edge > 2
% 1. fine boundary nodes >= 2 on this edge, phi as usual
[phi_1b, seg_num] = bezier_one_edge(nb_edge, nbridge_each);
ncbn_1b = seg_num * 2 + nbridge_each + 2;
elseif nb_edge == 2
% 2. fine boundary node = 1 on this edge, phi with 1 contro node
% linear interpolation
phi_1b = B2;
ncbn_1b = 3;
else
% 3. fine boundary node = 0 on this edge, phi = ones(2,2)
phi_1b = eye(4);
ncbn_1b = 2;
end
% to remove the adjoint corner_node
if edge == nedges
phi_1v(nbdofs : nbdofs + nb_edge*2 -1, ncbn : ncbn + (ncbn_1b-1)*2 - 1) = phi_1b(1:end-2, 1:end-2);
phi_1v(nbdofs : nbdofs + nb_edge*2 -1, 1:2) = phi_1b(1:end-2, end-1:end);
else
phi_1v(nbdofs : nbdofs + (nb_edge+1)*2 -1, ncbn : ncbn + ncbn_1b*2 - 1) = phi_1b;
end
% to remove the adjoint corner_node
ncbn = ncbn + (ncbn_1b-1) * 2;
nbdofs = nbdofs + nb_edge * 2;
end
phi_cell{ele} = phi_1v;
end
end
function B2 = simple_quad_Bezier(nodes_num)
B = zeros(2, 6, nodes_num);
for i = 0 : nodes_num-1
t = i / (nodes_num-1); % project to [0,1]
B(:, :, i+1) = Bezier_matrix(t);
end
B1 = permute(B, [2,1,3]); % [6, 2, x+1]
B2 = reshape(B1, 6, [])'; % [2(x+1), 6]
end
function B = Bezier_matrix(t)
b = [1, t, t^2] * [1, 0,0;
-2,2,0;
1,-2,1]; % 1,3
I = eye(2, 2);
B = kron(b, I); % 2,6
end