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.
84 lines
2.9 KiB
84 lines
2.9 KiB
3 years ago
|
%--------------------------
|
||
|
|
% @Author: Jingqiao Hu
|
||
|
|
% @Date: 2021-12-07 09:13:42
|
||
|
|
% @LastEditTime: 2022-01-13 21:12:16
|
||
|
|
|
||
|
|
% sensitivity on align-CVT of MMC
|
||
|
|
% dcdp = dcdH * dHdp
|
||
|
|
%--------------------------
|
||
|
|
function [c_dH, c_interEles] = sensitivity_mmc_align(nodes_cell, faces_cell, seeds0, step, T, ...
|
||
|
|
bpoly, epsilon, conn_list, c_rho, c_phi, c_interPolySeed, maxx, minx, maxy, miny)
|
||
|
|
|
||
|
|
npnts = size(seeds0, 1);
|
||
|
|
c_dH = cell(npnts, 3);
|
||
|
|
c_interEles = cell(npnts, 3);
|
||
|
|
|
||
|
|
seeds = [seeds0, T];
|
||
|
|
|
||
|
|
parfor i = 1 : npnts
|
||
|
|
|
||
|
|
neighbor = conn_list{i};
|
||
|
|
local_idx = int32([neighbor; i]);
|
||
|
|
|
||
|
|
for j = 1:3
|
||
|
|
% tic
|
||
|
|
% minus
|
||
|
|
seeds_minus = seeds;
|
||
|
|
seeds_minus(i, j) = seeds_minus(i, j) - step(j);
|
||
|
|
|
||
|
|
if j == 1 && seeds_minus(i, j) < minx
|
||
|
|
seeds_minus(i, j) = minx;
|
||
|
|
elseif j == 2 && seeds_minus(i, j) < miny
|
||
|
|
seeds_minus(i, j) = miny;
|
||
|
|
end
|
||
|
|
|
||
|
|
% only update the neighbor cells
|
||
|
|
[~, edges, c_nnbEdgeIDX, nnbSeedIDX] = intersection_info(seeds_minus(:,1:2), ...
|
||
|
|
local_idx, bpoly);
|
||
|
|
[MMCs, ~] = update_MMCs(edges, seeds_minus(:,3), nnbSeedIDX, seeds_minus(:,1:2));
|
||
|
|
|
||
|
|
[rhoMinus, interEles] = update_projection(nodes_cell, faces_cell, ...
|
||
|
|
MMCs, epsilon, c_rho, c_phi, c_nnbEdgeIDX, c_interPolySeed, local_idx);
|
||
|
|
|
||
|
|
% plus
|
||
|
|
seeds_plus = seeds;
|
||
|
|
seeds_plus(i, j) = seeds_plus(i, j) + step(j);
|
||
|
|
|
||
|
|
if j == 1 && seeds_plus(i, j) > maxx
|
||
|
|
seeds_plus(i, j) = maxx;
|
||
|
|
elseif j == 2 && seeds_plus(i, j) > maxy
|
||
|
|
seeds_plus(i, j) = maxy;
|
||
|
|
end
|
||
|
|
|
||
|
|
[~, edges, c_nnbEdgeIDX, nnbSeedIDX] = intersection_info(seeds_plus(:,1:2), ...
|
||
|
|
local_idx, bpoly);
|
||
|
|
[MMCs, ~] = update_MMCs(edges, seeds_plus(:,3), nnbSeedIDX, seeds_plus(:,1:2));
|
||
|
|
|
||
|
|
[rhoPlus, interEles] = update_projection(nodes_cell, faces_cell, ...
|
||
|
|
MMCs, epsilon, c_rho, c_phi, c_nnbEdgeIDX, c_interPolySeed, local_idx);
|
||
|
|
|
||
|
|
nnb = length(interEles);
|
||
|
|
dH_ele = cell(nnb, 1);
|
||
|
|
for k = 1 : nnb
|
||
|
|
ele = interEles(k);
|
||
|
|
dH_ele{k} = (rhoPlus{ele} - rhoMinus{ele}) / 2 / step(j);
|
||
|
|
% max(dH_ele{k}(:))
|
||
|
|
end
|
||
|
|
c_dH{i,j} = dH_ele;
|
||
|
|
c_interEles{i,j} = interEles;
|
||
|
|
|
||
|
|
% mp = [];
|
||
|
|
% figure;
|
||
|
|
% for testi = 1:size(MMCs_poly,1)
|
||
|
|
% mp = [mp; polyshape(reshape(MMCs_poly(testi,:), 2, [])')];
|
||
|
|
% end
|
||
|
|
% plot(mp);
|
||
|
|
% bp = [];
|
||
|
|
% for testi = 1:nnb
|
||
|
|
% bp = [bp; polyshape(c_polys{interEles(testi)})];
|
||
|
|
% end
|
||
|
|
% hold on; plot(bp);
|
||
|
|
% toc
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|