www.pudn.com > level_set_methods_1.1.zip > upwind_WENO.m, change:2005-05-05,size:2023b


function [data_x] = upwind_WENO(data, F, dx) 
% 
% Calculates the upwind derivative using 
% fifth order accurate WENO scheme 
% takes 1-D data 
% data: input data 
% F: force field 
% dx: grid resolution 
% Note: before entering this function, data and F need to be  
% extended by 3 at the beginning and end (values don't matter) 
% 
% Author: Baris Sumengen  sumengen@ece.ucsb.edu 
% http://vision.ece.ucsb.edu/~sumengen/ 
% 
 
 
if length(F) == 1 
    F = F*ones(size(data)); 
end 
 
if length(data) ~= length(F) 
    error('force and data lengths does not match'); 
end 
 
data_x = zeros(size(data)); 
 
% extrapolate the beginning and end points of data 
data(3) = 2*data(4)-data(5); 
data(2) = 2*data(3)-data(4); 
data(1) = 2*data(2)-data(3); 
data(end-2) = 2*data(end-3)-data(end-4); 
data(end-1) = 2*data(end-2)-data(end-3); 
data(end) = 2*data(end-1)-data(end-2); 
 
D1 = (data(2:end)-data(1:end-1))/dx; 
 
for i=1:(length(data)-6) 
    if F(i+3) > 0  % use D- 
        k = i-1; 
        v1 = D1(k+1); 
        v2 = D1(k+2); 
        v3 = D1(k+3); 
        v4 = D1(k+4); 
        v5 = D1(k+5); 
    elseif F(i+3) < 0  % use D+ 
        k = i; 
        v1 = D1(k+5); 
        v2 = D1(k+4); 
        v3 = D1(k+3); 
        v4 = D1(k+2); 
        v5 = D1(k+1); 
    else 
        continue; 
    end 
 
    data_x_1 = v1/3 - 7*v2/6 + 11*v3/6; 
    data_x_2 = -v2/6 + 5*v3/6 + v4/3; 
    data_x_3 = v3/3 + 5*v4/6 - v5/6; 
 
    epsilon = 1e-6*max([v1 v2 v3 v4 v5].^2) ... 
        + 1e-99; 
 
    S1 = (13/12)*(v1-2*v2+v3).^2 ... 
        + 0.25*(v1-4*v2+3*v3).^2; 
    S2 = (13/12)*(v2-2*v3+v4).^2 ... 
        + 0.25*(v2-v4).^2; 
    S3 = (13/12)*(v3-2*v4+v5).^2 ... 
        + 0.25*(3*v3-4*v4+v5).^2; 
 
    alpha1 = 0.1/((S1+epsilon).^2); 
    alpha2 = 0.6/((S2+epsilon).^2); 
    alpha3 = 0.3/((S3+epsilon).^2); 
    a_total = alpha1+alpha2+alpha3; 
    data_x(i+3) = (alpha1/a_total)*data_x_1 + (alpha2/a_total)*data_x_2 ... 
        + (alpha3/a_total)*data_x_3; 
 
    data_x(i+3) = F(i+3)*data_x(i+3); 
end