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.
56 lines
1.8 KiB
56 lines
1.8 KiB
function [H, diffH, phi_max] = projection_mmc(coorx, coory, MMCs, dx, cnum, H_filter, Hs_filter)
|
|
globaly = size(coorx, 1) - 1;
|
|
globalx = size(coorx, 2) - 1;
|
|
|
|
phi_max = repmat(-1e5, size(coorx));
|
|
phi_cell = cell(cnum, 1);
|
|
|
|
for ci = 1 : cnum
|
|
mmc = MMCs(:, ci);
|
|
phi = tPhi(mmc, coorx, coory);
|
|
phi_max = max(phi_max, phi);
|
|
phi_cell{ci} = phi;
|
|
end
|
|
|
|
% phi_max(:) = (H_filter*phi_max(:)) ./ Hs_filter;
|
|
|
|
epsilon = 4 * dx;
|
|
H = Heaviside(phi_max, globalx, globaly, epsilon);
|
|
H(:) = (H_filter*H(:)) ./ Hs_filter;
|
|
|
|
% sensitivity
|
|
step = [repmat(0.01, 1, 3), repmat(max(2*dx, 0.005), 1, 4)];
|
|
diffH = zeros(cnum, (globaly+1)*(globalx+1));
|
|
for ci = 1 : cnum
|
|
mmc = MMCs(:, ci);
|
|
|
|
mmc_plus = mmc;
|
|
mmc_plus(1) = mmc(1) + step(1);
|
|
tmpPhiD1 = tPhi(mmc_plus, coorx, coory);
|
|
phi_plus = tmpPhiD1;
|
|
for ik=1:ci-1
|
|
phi_plus = max(phi_plus, phi_cell{ik});
|
|
end
|
|
for ik=ci+1:cnum
|
|
phi_plus=max(phi_plus, phi_cell{ik});
|
|
end
|
|
|
|
mmc_minus = mmc;
|
|
mmc_minus(1) = mmc(1) - step(1);
|
|
tmpPhiD2 = tPhi(mmc_minus, coorx, coory);
|
|
phi_minus=tmpPhiD2;
|
|
for ik=1:ci-1
|
|
phi_minus=max(phi_minus, phi_cell{ik});
|
|
end
|
|
for ik=ci+1:cnum
|
|
phi_minus=max(phi_minus, phi_cell{ik});
|
|
end
|
|
% phi_plus(:) = (H_filter*phi_plus(:)) ./ Hs_filter;
|
|
% phi_minus(:) = (H_filter*phi_minus(:)) ./ Hs_filter;
|
|
HD1 = Heaviside(phi_plus, globalx,globaly,epsilon);
|
|
HD2 = Heaviside(phi_minus, globalx,globaly,epsilon);
|
|
HD1(:) = (H_filter*HD1(:)) ./ Hs_filter;
|
|
HD2(:) = (H_filter*HD2(:)) ./ Hs_filter;
|
|
diffH(ci, :) = (HD1-HD2) / (2*step(1));
|
|
end
|
|
end
|