A_operator.m

function [v,index_n]=A_operator(mu,u,index,Phi,M,epsilon,tau,dx)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                A-Operator fuer Buckley-Leverett Problem             %%%
%%%  A = Id-epsilon^2 tau D(u)_x d/dx - epsilon^2 tau D(u) d^2/dx^2     %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Input: mu        Parameter
%        u         Lsg u, wird komplett eingelesen und liest dann
%                  \in IR^{|index|x1} aus
%        index     Inidice-Menge, auf der lokale Lsg berechnet werden soll, 
%                  fuer index=[] global berechnen
%        Phi       Menge aus Basisvektoren, wird komplett eingelesen und 
%                  liest dann \in IR^{|index|xN} aus
%        M         Viskositaetsverhaeltnis
%        epsilon   Regularisierungsparameter
%        tau       Parameter
%        dx        Ortsschrittweite
% 
% Output: v        Auswertung des A-Operators, d.h. linke Seite des
%                  BL-Problems an Basisvektoren aus Phi, \in IR^{|index_n|xN}
%         index_n  neue Indexmenge, nimmt von jedem Miniintervall einen 
%                  Randpunkt weg
%


if isempty(index)  % globale Lsg
     n_x = length(u);
     a = zeros(n_x,1);
     c = zeros(n_x-1,1);
     for i=2:n_x
         a(i) = dx^2/(epsilon^2*tau)+BL_D(u(i-1),M)+BL_D(u(i),M);
     end
     a(1)=1;
     a(end)=1;
     for i=1:n_x-1
         c(i) = -BL_D(u(i),M);
     end
     A = sparse(epsilon^2*tau/dx^2*(diag(a)+diag(c,1)+diag(c,-1)));
     A(1,2);
     A(n_x,n_x-1);
     
     v = A*Phi;
    
     index_n = [];
    
else  % lokale Lsg
    
    index = sort(index);
    V = [];
    Ie = [];
    Ia = 1;
    I = [];
   
    for i=1:length(index)-1;
        if index(i+1)-index(i)==0
            error('Index kommt doppelt vor');
        elseif index(i+1)-index(i)>1     % Zusaetzliche Abfrage, dass mind 3 aufeinanderfolgende Werte
            Ia = [Ia, i+1];
            Ie = [Ie, i];
        end
    end
    Ie = [Ie,length(index)];
        
    for i=1:length(Ie)
        ni = Ie(i)-Ia(i)+1;  % Anzahl aufeinanderfolgender Stuetzstellen
        if ni<=2
           error('Zu wenig aufeinanderfolgende Stuetzstellen.') 
        end
        
        a = zeros(ni,1);
        c = zeros(ni-1,1);
        uc = u(index(Ia(i)):index(Ie(i)));
        for k=2:ni
            a(k) = dx^2/(epsilon^2*tau)+BL_D(uc(k-1),M)+BL_D(uc(k),M);
        end
        
        for k=1:ni-1
            c(k) = -BL_D(uc(k),M);
        end
        
        A = sparse(epsilon^2*tau/dx^2*(diag(a)+diag(c,1)+diag(c,-1)));
        vi = A*Phi(index(Ia(i)):index(Ie(i)),:);
        vi(1,:) = [];
        vi(end,:) = [];
        
        V = [V;vi];
        I = [I,index(Ia(i)+1):index(Ie(i)-1)];
        
    end
    v = V;
    index_n = I;
    
    
end