compute endlich hinzugefügt

kappa3 implementiert und in docu beschrieben

diff --git a/doc/doc.pdf b/doc/doc.pdf
index fa39c80304867e548f1e7d5adf19d489a8e051b9..49236f0cd75854352e51e92d38eec8b67934070c 100644 (file)
Binary files a/doc/doc.pdf and b/doc/doc.pdf differ

diff --git a/doc/doc.tex b/doc/doc.tex
index 21469dba3f907c8ce5e1359f7b8d028e3091a1d7..18da7854d847353744f138037136c65e073d3f95 100644 (file)
--- a/doc/doc.tex
+++ b/doc/doc.tex
@@ -658,9 +658,10 @@ Zum Plotten (\ref{exmplAA_2DQuad})werden noch folgende Schritte ausgeführt
  \item $error_{i} =  \sqrt{\enorm{\phi}^2 - \enorm{\phi_{l}}^2}$
  \item $\mu_{i} = \norm{\varrho^{1/2}(\phi_{\ell/2} - \phi_{l} )}$
  \item $\eta_{i} = \enorm{\phi_{\ell/2} - \phi_{l}}$
- \item $\kappa_{i} = \sqrt{\enorm{\phi_{\ell/2}}-\enorm{\phi_{\ell/2}^{(i-1)}}}$
- \item $\kappa2_{i} = \sqrt{\enorm{\phi_{l}}-\enorm{\phi_{l}^{(i-1)}}}$
-\end{itemize}
+ \item $\kappa_{i} = \enorm{\phi_{\ell/2}^{(i)}}-\enorm{\phi_{\ell/2}^{(i-1)}}$
+ \item $\kappa2_{i} = \enorm{\phi_{l}^{(i)}}-\enorm{\phi_{l}^{(i-1)}}$
+ \item $\kappa3_{i} = \enorm{\phi_{\ell/2}^{(i)}-\phi_{\ell/2}^{(i-1)}}$
+\end{itemize} 
 
 \begin{figure}[ht]
 \caption{2D Quad adaptiv anisotrop vollanalytisch $V\phi = 1$}

diff --git a/res/ebner2012.pdf b/res/ebner2012.pdf
new file mode 100755 (executable)
index 0000000..a8d2119
Binary files /dev/null and b/res/ebner2012.pdf differ

diff --git a/src/A_plots.m b/src/A_plots.m
index 9c4923ce0e61f4590b496d48e2ad3af23367c744..aa7e0d7fee3b854bacbdd9808798efaa8a4b54a1 100644 (file)
--- a/src/A_plots.m
+++ b/src/A_plots.m
@@ -12,7 +12,7 @@ leg2 = {};
 leg3 = {};
 sym = {};
 
-rows = 11;
+rows = 12;
 
 
 for i = 1:length(files)
@@ -52,6 +52,7 @@ for i = 1:length(files)
           ['\mu ' l0 l1{i}]...
           ['\kappa ' l0 l1{i}]...
           ['\kappa2 ' l0 l1{i}]...
+          ['\kappa3 ' l0 l1{i}]...
           }';
       leg1 = {leg1{:}...
           [ l0 l1{i}]...
@@ -100,23 +101,25 @@ i=0;
 % shift2 = shift2+shift2/10;
 % error*(eta(l)-shift2)/error(l)
 
-loglog(repmat(X(:,i+1),1,7),[G_D(:,2+rows*i) ...
+loglog(repmat(X(:,i+1),1,8),[G_D(:,2+rows*i) ...
     G_D(:,2+8+rows*i)...
     G_D(:,2+1+rows*i)...*(G_D(k,2)-shift)/G_D(k,3)...*G_D(1,2)/G_D(1,2+1+rows*i) ...
     sqrt(abs(sol - G_D(:,2+2+rows*i)))...*(G_D(k,2)-shift)/G_D(k,3)...
     G_D(:,2+3+rows*i)...*G_D(1,2)/G_D(1,2+3+rows*i) ...
     [ 0; sqrt(G_D(2:end,2+9+rows*i)-G_D(1:end-1,2+9+rows*i))]...
     [ 0; sqrt(G_D(2:end,2+2+rows*i)-G_D(1:end-1,2+2+rows*i))]...
+    G_D(:,2+10+rows*i)...
     ],type2sym{i+1});
 hold on
 for i = 1:step-1
-loglog(repmat(X(:,i+1),1,7),[G_D(:,2+rows*i) ...
+loglog(repmat(X(:,i+1),1,8),[G_D(:,2+rows*i) ...
     G_D(:,2+8+rows*i)...
     G_D(:,2+1+rows*i)...*G_D(1,2)/G_D(1,2+1+rows*i) ...
     sqrt(abs(sol - G_D(:,2+2+rows*i)))...*G_D(1,2)/sqrt(abs(sol - G_D(1,2+2+rows*i)))...
     G_D(:,2+3+rows*i)...*G_D(1,2)/G_D(1,2+3+rows*i) ...
     [ 0; sqrt(G_D(2:end,2+9+rows*i)-G_D(1:end-1,2+9+rows*i))]...
     [ 0; sqrt(G_D(2:end,2+2+rows*i)-G_D(1:end-1,2+2+rows*i))]...
+    G_D(:,2+10+rows*i)...
     ],type2sym{i+1});
     
 end

diff --git a/src/compute.m b/src/compute.m
new file mode 100644 (file)
index 0000000..8a79770
--- /dev/null
+++ b/src/compute.m
@@ -0,0 +1,274 @@
+function [data er] = compute(file,times,zeta,type,theta,nu,vcon,out)
+% [data er] = compute(file,times,zeta,type,theta,nu,vcon,out)
+% Führt times Verfeinerungsschritte aus
+%
+% file - StartNetz
+% times - wie oft Verfeinert werden soll
+% zeta & type - Bestimmen die Art der Matrix Berechnung
+% theta - adaptiv?
+% nu - isotrop?
+% vcon - Vorkonditionierung der Matrix? 1 oder 0
+% out - (optional) Dateizusatz um Netz, Lösung & co zu speichern
+%
+% P. Schaefer
+
+% Datei laden
+load(file)
+
+if(~exist('data','var'))
+    data=[];
+end
+kap3 = 0;
+%times mal verfeinern
+for j = 1:times  
+  if(exist('coo_fine','var'))
+    old_C_fine = coo_fine;
+    old_F_fine = f2s;
+    old_E_fine = ele_fine;
+    old_x_fine = x_fine;
+  end
+  
+  
+
+  %uniformIsotrop Verfeinern
+  [coo_fine,ele_fine,neigh_fine,f2s,sit_fine]=refineQuad(coordinates,elements,neigh,sites,2);
+
+  %Flaecheninhalte Berechnen (rhs)
+  b_fine = areaQuad(sit_fine);
+  
+  b = areaQuad(sites);
+  hmin = 2.^-max(sites,[],2);
+  hmax = 2.^-min(sites,[],2);
+
+
+  %data -> ErgebnisVektor
+  dataS = size(elements,1);
+  
+  %Alle MatrixBrechenungsArten mit dem selben Netz berechnen
+  for i = 1:length(type)
+    disp (size(elements,1))
+    %Matrix aufbauen -> MEX
+    V_fine = mex_build_V(coo_fine,ele_fine,zeta,type(i));
+    
+    %Testet auf Fehlerhafte Einträge (NaN +/-Inf)
+    [r c] = find(isnan(V_fine)~=isinf(V_fine));
+    if(~isempty(r))
+        figure(9)
+        plotShape(coo_fine,ele_fine(unique([r c]),:),'')
+        plotMark([r';c'],coo_fine,ele_fine)
+        title('Fehlerhafte Elemente')
+    end
+ 
+    if(~vcon)
+    %Lösung Berechnen
+        x_fine = V_fine\b_fine;
+        con = cond(V_fine);
+    else
+    %Vorkonditionierte Lösung!
+        D = diag(diag(V_fine.^(-1/2)));
+
+        A = D * V_fine * D;
+        c = D*b_fine;
+        y = A\c;
+        x_fine = D*y;
+        con = cond(A);
+    end
+    
+    % \tilde \mu ( \Pi h -h + L_2 )
+    tmu = hmin.* b .* sum((x_fine(f2s)'-repmat(sum(x_fine(f2s)',1)/4,4,1)).^2)' /4;
+    
+    %Fehlerschätzer 2 aufbauen
+    V = mex_build_V(coordinates,elements,zeta,type(i));
+    
+    if(~vcon)
+        x = V\b;
+    else
+        D = diag(diag(V.^(-1/2)));
+
+        A = D * V * D;
+        c = D*b;
+        y = A\c;
+        x = D*y;
+    end
+    
+    xo_fine(f2s) = repmat(x,1,4);
+    xd_fine = abs(xo_fine'-x_fine);
+    
+
+    
+    % \tilde \mu ( h/2 -h + L_2 )
+    mu = hmin.*b.*sum((x_fine(f2s)'-repmat(x',4,1)).^2)'/4;
+
+
+    %Energienorm^2 Berechnen |||h||| & |||h/2|||
+%     xe_fine = x_fine'*A_fine*x_fine;
+    xe_fine = b_fine'*x_fine;
+%     xe = x'*A*x;
+    xe = b'*x;
+    
+    
+    %Enorm^2 Elementweise Vergleich
+    
+    if(exist('old_F_fine','var'))
+%       xf_fine = zeros(size(x_fine));
+      for k = 1:size(ele_fine,1)
+        e_f = find(sum(k==f2s,2));
+        e_f_p = find(f2s(e_f,:)==k);
+        assert(size(e_f_p,1) ==1,'Error: e_f_p wurde mehrfach gefunden')
+        e_ff = find(sum(e_f==f,2));
+        e_ff_p = find(f(e_ff,:)==e_f);
+        assert(size(e_ff_p,1)==1,'Error: e_ff_p wurde mehrfach gefunden')
+        e_ff_s = length(e_ff_p);
+        if e_ff_s==4
+           e_p = e_f_p;
+        elseif e_ff_s == 1
+           e_p = e_ff_p;  
+        elseif e_ff_s == 2      
+           if(sum(e_ff_p) == 5)
+             e_p = e_ff_p(floor((e_f_p-1)/2)+1);
+           else
+             inh = [1 2 4 3];
+             e_p = inh(e_ff_p(mod(floor((e_f_p)/2),2)+1));
+           end
+        end
+%         figure(11)
+%           plotShape(old_C,old_E,'b');view(2);
+%           plotMark(e_ff,coo_fine,ele_fine);     
+%           
+%         figure(12)
+%           plotShape(old_C_fine,old_E_fine,'b'); view(2);
+%           plotMark(old_F_fine(e_ff,e_p),coo_fine,ele_fine);
+%           
+%         
+%         figure(13)
+%           plotShape(coordinates,elements,'b');view(2);
+%           plotMark(e_f,coo_fine,ele_fine);
+% 
+%         
+%         figure(14)
+%           plotShape(coo_fine,ele_fine,'b');view(2);
+%           plotMark(k,coo_fine,ele_fine);
+        
+        xf_fine(k,1) = abs(x_fine(k) - old_x_fine(old_F_fine(e_ff,e_p)));
+      end
+
+      kap3 = b_fine'*(xf_fine);
+    end
+      
+%     xe_fine = b_fine'*x_fine;
+
+
+    
+    %\tilde \mu 2 = ( ||\Pi h|| - ||h||)
+    tmu2 = hmin.* b.* (...
+    sum((x_fine(f2s)').^2)'-sum(repmat(sum(x_fine(f2s)',1)/4,4,1).^2)'...
+    ) /4;
+  
+    % |||h/2 -h|||
+%     eta = xd_fine'*A_fine*xd_fine;    
+    eta = abs(xe_fine-xe);
+    
+%    save_A_fine{i} = V_fine;
+%    save_x_fine{i} = x_fine;
+%    
+%    save_A{i} = V;
+%    save_x{i} = x;
+   
+    
+    
+    
+    dataS = [dataS ...
+        type(i) ... berechnet mit Typ
+        sqrt(sum(tmu))... tilde mu
+        sqrt(eta) ... eta
+        xe ... error (kappa 2)
+        sqrt(sum(mu))... mu
+        min(hmin)/max(hmax)...
+        min(hmax)/max(hmax)...
+        min(hmin./hmax) con...
+        sqrt(sum(tmu2))... tilde mu 2
+        xe_fine... (kappa)
+        kap3 ... kappa 3
+        ];
+  end
+
+  
+  % nur RandElemente Verfeinern
+%   marked = ones(1,size(G_E,1));
+%   marked(find(sum((G_N(:,1:4)==0),2))) = 2;
+
+  % Markieren mit gewählten Parametern
+  marked = mark(x_fine(f2s)',tmu,theta,nu);
+  
+  % Netz bunt Plotten!
+%   figure(1)
+%   plotShape(G_C,G_E,'s',tmu);
+%   title('Elemente mit Fehlerschaetzer')
+%   colorbar
+%   view(2)
+%   plotMark(find(marked>1),G_C,G_E);
+  
+  assert(size(elements,1)==length(marked),'MarkierungsVektor ist fehlerhaft')
+  
+  old_C = coordinates;
+  old_E = elements;
+  old_S = sites;
+  
+%   figure(1)
+%   plotShape(coordinates,elements(:,:),'db');view(2);
+    
+  %Netz Verfeinern, wie durch marked bestimmt
+  [coordinates, elements, neigh, f, sites, er] = refineQuad(coordinates,elements,neigh,sites,marked);
+  
+%   figure(2)
+%   plotShape(coordinates,elements(:,:),'db');view(2);
+  
+%   f
+  
+  %Vater Sohn test
+  assert(sum(areaQuad(old_S))==sum(areaQuad(sites)),'Gesamtinhalt Fehlerhaft')
+  
+  e = 0;
+  for k = 1:size(data,1)
+    e = e +  abs(areaQuad(old_S(k,:)) - sum(areaQuad(sites(unique(f(k,:)'),:))));
+  end
+  assert(e==0, 'Einzelne Inhalte Fehlerhaft')
+  
+  % Fehlerhafte Elemente Plotten
+%   if(~isempty(er))
+%     figure(10)
+%     plotShape(G_C,G_E,'');
+%     plotMark(er,G_C,G_E);
+%   end
+ 
+  %ErgebnisWerte Speichern
+  data(size(data,1)+1,1:length(dataS)) = dataS;
+  typeN = int2str(type);
+  save (['meshSave/' out typeN(typeN~=' ') '_' int2str(size(data,1))], 'coordinates', 'elements','neigh','data', 'sites')
+  
+  %Alle Relevanten zwischenInformationen Speichern
+%   out = '_';
+%   if(length(varargin)~=0)
+%      out = varargin{1};
+%    end
+%   typeN = int2str(type);
+%  save (['meshSave/fine_' out typeN(typeN~=' ') '_' int2str(size(data,1))],...
+%      'coo_fine', 'ele_fine','neigh_fine','f2s','data',...
+%      'save_A','save_x','save_A_fine','save_x_fine','b','b_fine')
+
+%   clear 'coo_fine' 'ele_fine' 'neigh_fine' 'f2s' 
+%   plotShape(G_C,G_E,'');
+
+end
+end
+
+
+% kappa:
+% kappa = ||| phi_{h/2}^{\ell} - \phi_{h/2}^{\ell-1} |||
+% ist die differenz der eNormen von den h/2 Netzen?
+% 
+% tildeMu:
+% tilde-mu(T)^2 = hmin(T) || (1-Pi_h) phi_{h/2} ||^2_{L2(T)}
+% = hmin(T) * ( ||phi_{h/2}||^2_{L2(T)} - ||Pi_h
+% phi_{h/2}||^2_{L2(T)}
+