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
% 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);
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
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;
% 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')
% 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)
% 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
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)';
% 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)
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));
%
% 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
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;
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;
%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;
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);
% 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
end
+ % Da Nachbarn noch Verfeinert werden müssen erst mal weiter
continue;
end
end
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
%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;
%% 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)],:);
%% 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;
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);
% %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
projects / bacc.git / commitdiff