From 044151ef6a103e1d47ac28649adb1bc82348be32 Mon Sep 17 00:00:00 2001 From: treecity Date: Fri, 2 Mar 2012 21:30:29 +0000 Subject: [PATCH] [src] Wichtige Funktionen kommentiert\ [src] mark Isotrop wird nun richtig abgefangen git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@93 26120e32-c555-405d-b3e1-1f783fb42516 --- src/A_plots.m | 7 ++-- src/A_run.m | 27 ++++++++++------ src/A_step.m | 65 ++++++++++++++++++++++++++------------ src/computeEstSlpMuTilde.m | 6 +++- src/mark.m | 10 +++--- src/refineQuad.m | 50 ++++++++++++++++++----------- 6 files changed, 107 insertions(+), 58 deletions(-) diff --git a/src/A_plots.m b/src/A_plots.m index 5e45a34..c77e4c5 100644 --- a/src/A_plots.m +++ b/src/A_plots.m @@ -56,11 +56,12 @@ end if step<1 disp ('Error: No Data to show.') else - G_D + format longG + data((end-8):end,[1 3 4]) - sol = interp1(1./X((end-5):(end-1),3)',G_D((end-5):(end-1),3)',0,'spline') + sol = interp1(1./X((1):(end-2),3)',G_D((1):(end-2),3)',0,'spline') % sol = 8.28466; - + format short figure(4) loglog(X(:,[2+(0:step-1)*3]),G_D(:,[2+(0:step-1)*3]),'--o') hold on diff --git a/src/A_run.m b/src/A_run.m index 216084e..d3a5131 100644 --- a/src/A_run.m +++ b/src/A_run.m @@ -7,17 +7,24 @@ function A_run(file,times,mu,type,eta,eps,out) % eta - alle verfeinern oder nur wichtige? % eps - wie verfeinern iso or aniso % out - Dateiname der Ausgabe +% +% P. Schaefer +%Parameter testen +assert(length(mu)==max(type)-1||length(mu)==1,... + 'Pleas set right type and mu parameters'); +%Netz laden A_loadMeshF(file) +%globale Speicherung initialisieren global G_T; global G_D; -assert(length(mu)==max(type)-1||length(mu)==1,... - 'Pleas set right type and mu parameters'); - +%Schrittweise Verfeinern for i = 1:times + + %benötigte Zeit Schätzen if(size(G_T,1)>2) nextF = G_T(end,1)*4; nextS = neville2(1:size(G_T,1),G_T(:,1)',size(G_T,1)+1); @@ -25,16 +32,18 @@ for i = 1:times calc = size(G_T,1)+(-2:0); nextTime = [nextS neville2(G_T(calc,1)',G_T(calc,2)',nextF) ... neville2(G_T(cald,1)',G_T(cald,3)',nextF) ... - neville2(G_T(calc,1)',G_T(calc,4)',nextS)] - - %overFlowTime = neville2(G_T(size(G_T,1)+(-1:0),1)',sum(G_T(size(G_T,1)+(-1:0),2:4),2)',4500) + neville2(G_T(calc,1)',G_T(calc,4)',nextS)]; end - + + % Ein kompletter Verfeinerungschritt A_step(mu,type,eta,eps,out); - usedTime = G_T(end,:) + + %Zeit Speichern + usedTime = G_T(end,:); + + %Ergebnisse Abspeichern data = G_D(end,:) typeN = int2str(type); - A_saveMesh([out typeN(typeN~=' ') '_' int2str(size(G_T,1))]); end diff --git a/src/A_step.m b/src/A_step.m index d59cf53..df1de02 100644 --- a/src/A_step.m +++ b/src/A_step.m @@ -1,60 +1,76 @@ function [data time er] = A_step(mu,type,eta,eps,varargin) -% [dataAniso time er] = A_stepAnIso(mu,type,eta,eps) - +% [data time er] = A_step(mu,type,eta,eps [, out]) +% Führt einen Verfeinerungsschritt aus +% +% mu & type - Bestimmen die Art der Matrix Berechnung +% eta - adaptiv? +% eps - isotrop? +% out - (optional) Dateizusatz um Netz, Lösung & co zu speichern +% +% P. Schaefer + +%globale Speicherung initialisieren global G_E; global G_C; global G_N; global G_T; global G_D; +%Netz darf nicht leer sein assert(~isempty(G_E) && ~isempty(G_C) && ~isempty(G_N),... 'Please use A_loadMesh first') +%Steuerungs Parameter prüfen assert(length(mu)==max(type)-1||length(mu)==1,... 'Pleas set right type and mu parameters') +%Steuerungsparameter korrigieren if(length(mu)==1) mu = repmat(mu,max(type)-1,1); end -out = ''; -if(length(varargin)~=0) - out = varargin{1}; -end - +%Vektor für Zeitmessung anlegen time = zeros(1,3); tic + %uniformIsotrop Verfeinern [coordinates_fine,elements_fine,neigh_fine,f2s]=refineQuad(G_C,G_E,G_N,2); time(1) = toc; + %Flaecheninhalte Berechnen (rhs) b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2)); tic + %data -> ErgebnisMatrix data = size(G_E,1); + %save_* -> ZwischenSpeicherung save_A = {}; save_x = {}; + + %Alle MatrixBrechenungsArten mit dem selben Netz berechnen for i = 1:length(type) + %Matrix aufbauen -> MEX A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type(i)); - [r c] = find(isnan(A_fine)~=isinf(A_fine)); - if(~isempty(r)) + %Testet auf Fehlerhafte Einträge (NaN +/-Inf) + [r c] = find(isnan(A_fine)~=isinf(A_fine)); + if(~isempty(r)) figure(9) plotShape(coordinates_fine,elements_fine(unique([r c]),:),'') plotMark([r';c'],coordinates_fine,elements_fine) title('Fehlerhafte Elemente') - end + end - - + %Lösung Berechnen x_fine = A_fine\b; + %Energienorm^2 Berechnen xe_fine = x_fine'*A_fine*x_fine; + %Fehlerschätzer aufbauen ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s); save_A{i} = A_fine; save_x{i} = x_fine; - data = [data type(i) sqrt(sum(ind)) xe_fine]; end time(2) = toc; @@ -63,8 +79,10 @@ time = zeros(1,3); % marked = ones(1,size(G_E,1)); % marked(find(sum((G_N(:,1:4)==0),2))) = 2; - % Adaptiv Verfeinern + % Markieren mit gewählten Parametern marked = mark(x_fine(f2s)',ind,eta,eps); + + % Bunt Plotten! % figure(1) % plotShape(G_C,G_E,'s',ind); % title('Elemente mit Fehlerschaetzer') @@ -72,14 +90,14 @@ time = zeros(1,3); % view(2) % plotMark(find(marked>1),G_C,G_E); - -% clear 'coordinates_fine' 'elements_fine' 'neigh_fine' 'f2s' - assert(size(G_E,1)==length(marked),'MarkierungsVektor ist fehlerhaft') tic er = []; + + %Netz Verfeinern, wie durch marked bestimmt [G_C, G_E, G_N, f, er] = refineQuad(G_C,G_E,G_N,marked); + time(3) = toc; % if(~isempty(er)) % figure(10) @@ -87,16 +105,21 @@ time = zeros(1,3); % plotMark(er,G_C,G_E); % end - + %ErgebnisWerte und Zeiten Speichern G_T(size(G_T,1)+1,1:4) = [size(G_E,1) time]; G_D(size(G_D,1)+1,1:length(data)) = data; -% A_saveMesh([out typeN(typeN~=' ') '_' int2str(size(G_T,1))]); + + %Alle Relevanten zwischenInformationen Speichern + out = ''; + if(length(varargin)~=0) + out = varargin{1}; + end typeN = int2str(type); save (['meshSave/fine' out typeN(typeN~=' ') '_' int2str(size(G_T,1))],... 'coordinates_fine', 'elements_fine','neigh_fine','f2s','time','data','save_A','save_x') - - + +% clear 'coordinates_fine' 'elements_fine' 'neigh_fine' 'f2s' % plotShape(G_C,G_E,''); end diff --git a/src/computeEstSlpMuTilde.m b/src/computeEstSlpMuTilde.m index 78e40ce..86e481f 100644 --- a/src/computeEstSlpMuTilde.m +++ b/src/computeEstSlpMuTilde.m @@ -1,9 +1,13 @@ function ind=computeEstSlpMuTilde(x_fine,coo,ele,f2s) +% ind=computeEstSlpMuTilde(x_fine,coo,ele,f2s) +% Berechnet den Fehlerschätzer über die Elemente +% % x_fine - Phi % coo - coordinaten % ele - elemente % f2s - Vater Sohn beziehungen - +% +% P. Schaefer xF2S = x_fine(f2s)'; diff --git a/src/mark.m b/src/mark.m index 1d62366..61f5290 100644 --- a/src/mark.m +++ b/src/mark.m @@ -2,8 +2,8 @@ function REF = mark(xF2S,ind,theta,eta) % function REF = mark(xF2S,ind,theta,eta) % xF2S - Father son relation % ind - error estimator -% eta - refine element? (0..1, 1 = All) -% eps - refine how? (0 = Isotrop) +% theta - refine element? (0..1, 1 = All) +% eta - refine how? (0...1, 0 = Isotrop) % REF - vector with entries [1 : 4] if(size(xF2S,1)==1) @@ -26,15 +26,15 @@ if(theta <1) ell = find(sum_ind >= sum_ind(end) * theta,1); - %Symertrisieren + %Symmetrisieren ell = ell + find(abs(( sum_ind(ell)-sum_ind(ell:end)))/sum_ind(ell)>10^-2,1); - t1(idx(ell+1:end)) = 1; % Nicht verfeinert + t1(idx(ell+1:end)) = 1; % Nicht verfeinern end %% Wie muss verfeinert werden -if(eps > 0) % Horizontal oder Vertikal +if(eta > 0) % Horizontal oder Vertikal t3 = (eta*abs(Ct(3,:)) >= sqrt(Ct(2,:).^2+Ct(4,:).^2)); REF(t3-t1==1) = 3; t4 = (eta*abs(Ct(4,:)) >= sqrt(Ct(2,:).^2+Ct(3,:).^2)); diff --git a/src/refineQuad.m b/src/refineQuad.m index c6fe523..e784f10 100644 --- a/src/refineQuad.m +++ b/src/refineQuad.m @@ -4,13 +4,14 @@ function [coo,ele,nei,f2s,err] = refineQuad(coordinates,elements,neigh,type) % % Verfeinert die markierten Elemente mit dem entsprechenden TYP und gibt % auch die F2S Beziehungen zurück. type muss von der Länge der Anzahl der -% Elemente entsprechen oder genau 1 und die Einträge können 1,2,3,4 sein. +% Elemente entsprechen oder genau 1 und die Einträge können 1,2,3,4,5 sein. % % der Typ zu jedem Element entspricht dabei: -% 1- keine Verfeinerung -% 2- 4 neue Elemente -% 3- 2 neue Elemente, übereinander -% 4- 2 neue Elemente, nebeneinander +% 1 - keine Verfeinerung +% 2 - 4 neue Elemente +% 3 - 2 neue Elemente, übereinander +% 4 - 2 neue Elemente, nebeneinander +% 5 - 4 neue Elemente, wird erst (3) und dann beide Elemente (4) geteilt % % P. Schaefer @@ -24,6 +25,7 @@ end assert(size(elements,1)==size(neigh,1)&&size(elements,1)==length(type),... 'Dimensionen passen nicht'); + %Globale Variabelen aufbauen global G_ref_E; global G_ref_C; @@ -33,7 +35,7 @@ global G_ref_t; %wie soll verfeinert werden global G_ref_tD; %wie wurde bereits verfeinert (in diesem Durchlauf) global G_ref_f2sT; %Temporäre Beziehung -%Globale Variablen zuweisen +%INTERNE Globale Variablen zuweisen G_ref_E = elements; G_ref_C = coordinates; G_ref_N = neigh; @@ -55,11 +57,15 @@ ref_old2 = []; %Jedes Element verfeinern while(1==1) - %TODO jeden 4ten 2er durch schrittweise Verfeinerung ersetzen + %jeden 4ten 2er durch schrittweise Verfeinerung (5) ersetzen t_ref=find(G_ref_t==2); G_ref_t(t_ref(1:4:end)) = 5; + + %Welche Elemente müssen bearbeitet werden ref = find(G_ref_t>1); ref = reshape(ref,1,length(ref)); + + %Muss noch weiter Verfeinert werden? if(isequal(ref,ref_old)) assert(~isequal(ref_old2,G_ref_t,'Markierte sind verschieden')); break; @@ -69,16 +75,22 @@ while(1==1) if(isempty(ref)) break; end + + % figure(6) % plotShape(G_ref_C,G_ref_E) % plotMark(ref,G_ref_C,G_ref_E,'xg') % title('Zum Verfeinern Markierte Elemente') - for ele = ref % ref(randperm(length(ref))) %LOL Zufall bringts + + % Elementeweise Bearbeiten + for ele = ref % ref(randperm(length(ref))) + % figure(5) % plotShape(G_ref_C,G_ref_E) % plotMark(ele,G_ref_C,G_ref_E,'xg') % view(2) % title('Entscheidene Nachbarelemente') + % # HangingNode Check Nt = find(G_ref_N(ele,5:8)==0); N = G_ref_N(ele,Nt); @@ -95,6 +107,7 @@ while(1==1) % plotMark(N2,G_ref_C,G_ref_E,'or') + %Hat noch zu teilende Nachbarn? if(~isempty(N2)) N3t = mod(find((G_ref_N(N2',:)==ele)')-1,4)+1; %Nachbarseiten N4t = find(diag(G_ref_N(N2',(N3t + 4)'))~=0)'; %Nachbarn mit 2Nachbarn @@ -121,6 +134,7 @@ while(1==1) end + % Da Nachbarn noch Verfeinert werden müssen erst mal weiter continue; end end @@ -129,9 +143,9 @@ while(1==1) assert(G_ref_tD(ele)~=2,'Element ist schon verfeinert') assert(G_ref_t(ele)>1,'Element ist nicht Markiert') G_ref_f2sT = ones(1,4)*ele; - refineE(ele); - updateN(ele); - updateF2S(ele); + refineE(ele); %Element Teilen + updateN(ele); %Nachbarn des Elements aktualisieren + updateF2S(ele); %VaterSohn Relation setzen % plotShape(G_ref_C,G_ref_E) end end @@ -144,16 +158,16 @@ f2s = G_ref_f2s; %Doppelte Koordinaten loeschen [coo l pos] = unique(coo,'rows'); - ele = pos(ele); - %Globale Variablen freigeben + %INTERNE Globale Variablen freigeben clear G_ref_E G_ref_C G_ref_N G_ref_f2s G_ref_t G_ref_tD G_ref_s end %% Element Verfeinern ! sollte nur ausgeführt werden wenn wirklich möglich function refineE(ele) % Element wird gnadenlos Verfeinert + global G_ref_E; global G_ref_C; global G_ref_f2sT; @@ -202,13 +216,14 @@ end %% Aktualisieren der Nachbarn ! sollte nur ausgeführt werden wenn wirklich möglich function updateN(ele) % Nachbarschaften werden neu gesetzt (nach N und f2s) + global G_ref_E; global G_ref_N; global G_ref_f2sT; global G_ref_t; this = G_ref_N(ele,:); - %An welchen Kanten habe ich Nachbarn +%An welchen Kanten habe ich Nachbarn S = find(mod((this(1:4)~=0).*(this(5:8)==0),2))'; %Einen Nachbar (Single) D = find(this(5:8)~=0)'; %Zwei Nachbarn (Double) % G_ref_N([this(S) this(D)],:); @@ -321,6 +336,7 @@ end %% Aktualisieren der VaterSohn Beziehung ! sollte nur ausgeführt werden wenn wirklich möglich function updateF2S(ele) +%Vater Sohn Beziehungen richtig Setzen global G_ref_f2s; global G_ref_f2sT; global G_ref_t; @@ -334,13 +350,12 @@ function updateF2S(ele) if(G_ref_t(ele)<5) G_ref_tD(G_ref_f2sT) = G_ref_t(ele); G_ref_t(G_ref_f2sT) = 0; - else + else %Element muss noch mal geteilt werden G_ref_tD(G_ref_f2sT) = 3; G_ref_t(G_ref_f2sT) = 4; end else %Wenn Element zum zweiten Mal verfeinert wird G_ref_tD(G_ref_f2sT) = 2; - %TODO f2s AKTUALISIEREN org=floor((find(G_ref_f2s'==ele,1)-1)/4)+1; pos=find(G_ref_f2s(org,:)==ele,1); @@ -375,9 +390,6 @@ end % %Globale Variabelen aufbauen % global G_ref_E; % global G_ref_C; -% % global G_ref_N; -% % global G_ref_f2s; -% % global G_ref_t; % global G_ref_s; % for i = ele % % normalized Vector on every triangle -- 2.47.3