rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
b_operator.m
1 function [b,index_n]=b_operator(mu,u,index,M,epsilon,dx)
2 %
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 %%% b-Operator fuer Buckley-Leverett Problem %%%
5 %%% b = epsilon D(u)_x u_x + D(u)u_xx - f(u)_x %%%
6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 %
8 % Input: mu Parameter
9 % u Lsg u, wird komplett eingelesen und liest dann
10 % \in IR^{|index|x1} aus
11 % index Inidice-Menge, auf der lokale Lsg berechnet werden soll,
12 % fuer index=[] global berechnen
13 % M Viskositaetsverhaeltnis
14 % epsilon Regularisierungsparameter
15 % dx Ortsschrittweite
16 %
17 % Output: b Auswertung des b-Operators, d.h. rechte Seite des
18 % BL-Problems, \in IR^{|index_n|x1}
19 % index_n neue Indexmenge, nimmt von jedem Miniintervall einen
20 % Randpunkt weg
21 
22 %
23 
24 if isempty(index) % globale Lsg
25 
26  n_x = length(u);
27  b = zeros(n_x,1); % initialisieren
28  b(1) = epsilon*BL_D(u(1),M)*(u(2)-u(1))/(dx)^2; % OE Intervall so groß, inflow/outflow
29  for i=2:n_x-1
30  b(i) = epsilon*(BL_D(u(i),M)*(u(i+1)-u(i))-BL_D(u(i-1),M)*(u(i)-u(i-1)))/(dx^2)-(BL_f(u(i),M)-BL_f(u(i-1),M))/dx;
31  end
32  b(n_x) = -epsilon*BL_D(u(n_x-1),M)*(u(n_x)-u(n_x-1))/(dx)^2-(BL_f(u(n_x),M)-BL_f(u(n_x-1),M))/dx;
33  index_n = [];
34 
35 else % lokale Lsg
36  index = sort(index);
37  B = [];
38  Ie = [];
39  Ia = 1;
40  I = [];
41  for i=1:length(index)-1;
42  if index(i+1)-index(i)==0
43  error('Index kommt doppelt vor');
44  elseif index(i+1)-index(i)>1 % Zusaetzliche Abfrage, dass mind 3 aufeinanderfolgende Werte
45  Ia = [Ia, i+1];
46  Ie = [Ie, i];
47  end
48  end
49  Ie = [Ie,length(index)];
50 
51  for i=1:length(Ie)
52  ni = Ie(i)-Ia(i)+1; % Anzahl aufeinanderfolgender Stuetzstellen
53  if ni<=2
54  error('Zu wenig aufeinanderfolgende Stuetzstellen.')
55  end
56 
57  uc = u(index(Ia(i)):index(Ie(i)));
58  b = zeros(ni,1); % initialisieren
59  for k=2:ni-1
60  b(k) = epsilon*(BL_D(uc(k),M)*(uc(k+1)-uc(k))-BL_D(uc(k-1),M)*(uc(k)-uc(k-1)))/(dx^2)-(BL_f(uc(k),M)-BL_f(uc(k-1),M))/dx;
61  end
62  b(1) = [];
63  b(end) = [];
64 
65  B = [B;b];
66  I = [I,index(Ia(i)+1):index(Ie(i)-1)];
67 
68  end
69  b = B;
70  index_n = I;
71 
72 end