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.
102 lines
3.4 KiB
102 lines
3.4 KiB
3 years ago
|
%--------------------------
|
||
|
|
% @Author: Jingqiao Hu
|
||
|
|
% @Date: 2021-12-09 11:54:24
|
||
|
|
% @LastEditTime: 2021-12-15 11:00:41
|
||
|
|
|
||
|
|
% pb for blending
|
||
|
|
%--------------------------
|
||
|
|
function [rho_cell, dH_cell] = projection_mmc_cvt(nodes_cell, faces_cell, MMCs, epsilon, nd, pb)
|
||
|
|
|
||
|
|
% nd = 3;
|
||
|
|
|
||
|
|
ncvt = size(nodes_cell, 1);
|
||
|
|
|
||
|
|
rho_cell = cell(ncvt, 1);
|
||
|
|
dH_cell = cell(ncvt, 1);
|
||
|
|
|
||
|
|
% figure; clf
|
||
|
|
for ele = 1 : ncvt
|
||
|
|
nodes = nodes_cell{ele};
|
||
|
|
faces = faces_cell{ele};
|
||
|
|
|
||
|
|
mmc = MMCs{ele};
|
||
|
|
nc = size(mmc, 1);
|
||
|
|
|
||
|
|
phi_max = repmat(-1e5, size(nodes, 1), 1);
|
||
|
|
% phi0 = 0;
|
||
|
|
phi_cell = cell(nc, 1);
|
||
|
|
for i = 1 : nc
|
||
|
|
phi = tPhi(mmc(i, :), nodes(:,1), nodes(:,2));
|
||
|
|
phi_max = max(phi_max, phi);
|
||
|
|
phi_cell{i} = phi;
|
||
|
|
% phi0 = phi0 + exp(pb * phi);
|
||
|
|
end
|
||
|
|
% phi_max = log(phi0) / pb;
|
||
|
|
|
||
|
|
dx = 4 * epsilon;
|
||
|
|
H = Heaviside(phi_max, dx);
|
||
|
|
|
||
|
|
rho_cell{ele} = sum(H(faces), 2) / 4;
|
||
|
|
|
||
|
|
% patch('Faces',faces,'Vertices',nodes,'FaceVertexCData',1 - H','FaceColor','interp','EdgeColor','interp'); hold on;
|
||
|
|
|
||
|
|
%% sensitivity
|
||
|
|
step = repmat(epsilon / 2, 1, nd);
|
||
|
|
dH = zeros(nc, nd, size(faces, 1));
|
||
|
|
for i = 1 : nc
|
||
|
|
for di = 1 : nd
|
||
|
|
mmc_plus = mmc(i, :);
|
||
|
|
mmc_plus(di) = mmc(di) + step(di);
|
||
|
|
tmpPhiD1 = tPhi(mmc_plus, nodes(:,1), nodes(:,2));
|
||
|
|
phi_plus = tmpPhiD1;
|
||
|
|
% phi_plus = exp(pb * tmpPhiD1);
|
||
|
|
for ik = 1:i-1
|
||
|
|
% phi_plus = phi_plus + exp(pb * phi_cell{ik});
|
||
|
|
phi_plus = max(phi_plus, phi_cell{ik});
|
||
|
|
end
|
||
|
|
for ik = i+1:nc
|
||
|
|
% phi_plus = phi_plus + exp(pb * phi_cell{ik});
|
||
|
|
phi_plus = max(phi_plus, phi_cell{ik});
|
||
|
|
end
|
||
|
|
% phi_plus_blend = log(phi_plus) / pb;
|
||
|
|
|
||
|
|
mmc_minus = mmc(i, :);
|
||
|
|
mmc_minus(di) = mmc(di) - step(di);
|
||
|
|
tmpPhiD2 = tPhi(mmc_minus, nodes(:,1), nodes(:,2));
|
||
|
|
phi_minus=tmpPhiD2;
|
||
|
|
% phi_minus = exp(pb * tmpPhiD2);
|
||
|
|
for ik = 1:i-1
|
||
|
|
% phi_minus = phi_minus + exp(pb * phi_cell{ik});
|
||
|
|
phi_minus = max(phi_minus, phi_cell{ik});
|
||
|
|
end
|
||
|
|
for ik = i+1:nc
|
||
|
|
% phi_minus = phi_minus + exp(pb * phi_cell{ik});
|
||
|
|
phi_minus = max(phi_minus, phi_cell{ik});
|
||
|
|
end
|
||
|
|
% phi_minus_blend = log(phi_minus) / pb;
|
||
|
|
|
||
|
|
% HD1 = Heaviside(phi_plus_blend, dx);
|
||
|
|
% HD2 = Heaviside(phi_minus_blend, dx);
|
||
|
|
HD1 = Heaviside(phi_plus, dx);
|
||
|
|
HD2 = Heaviside(phi_minus, dx);
|
||
|
|
rho1 = sum(HD1(faces), 2) / 4;
|
||
|
|
rho2 = sum(HD2(faces), 2) / 4;
|
||
|
|
dH(i, di, :) = (rho1 - rho2) / (2 * step(di));
|
||
|
|
end
|
||
|
|
end
|
||
|
|
dH_cell{ele} = dH;
|
||
|
|
end
|
||
|
|
% colormap("gray"); axis equal; axis off; pause(1e-6);
|
||
|
|
end
|
||
|
|
|
||
|
|
function H = Heaviside(phi, epsilon)
|
||
|
|
alpha=1e-3; % parameter alpha in the Heaviside function
|
||
|
|
|
||
|
|
num_all=[1 : length(phi)]';
|
||
|
|
num1=find(phi>epsilon);
|
||
|
|
H(num1)=1;
|
||
|
|
num2=find(phi<-epsilon);
|
||
|
|
H(num2)=alpha;
|
||
|
|
num3=setdiff(num_all,[num1;num2]);
|
||
|
|
H(num3)=3*(1-alpha)/4*(phi(num3)/epsilon-phi(num3).^3/(3*(epsilon)^3))+ (1+alpha)/2;
|
||
|
|
end
|