From: Peter Schaefer Date: Thu, 19 Jul 2012 13:29:08 +0000 (+0200) Subject: compute endlich hinzugefügt X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=edeebfa5a12d3fc7ef1b39e0894ba3425c9e10e0;p=bacc.git compute endlich hinzugefügt kappa3 implementiert und in docu beschrieben --- diff --git a/doc/doc.pdf b/doc/doc.pdf index fa39c80..49236f0 100644 Binary files a/doc/doc.pdf and b/doc/doc.pdf differ diff --git a/doc/doc.tex b/doc/doc.tex index 21469db..18da785 100644 --- 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 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 9c4923c..aa7e0d7 100644 --- 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 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)} +