--- /dev/null
+function A_AnIso(file,times,mu,type,eta,eps)
+%A_ANISO Summary of this function goes here
+% Detailed explanation goes here
+
+A_loadMeshF(file)
+
+global G_Ta;
+global G_Da;
+
+for i = 1:times
+ if(size(G_Ta,1)>2)
+ nextF = G_Ta(end,1)*4;
+ nextS = neville2(1:size(G_Ta,1),G_Ta(:,1)',size(G_Ta,1)+1);
+ calc = size(G_Ta,1)+(-1:0);
+ nextTime = [nextS neville2(G_Ta(calc,1)',G_Ta(calc,2)',nextF) ...
+ neville2(G_Ta(calc,1)',G_Ta(calc,3)',nextF) ...
+ neville2(G_Ta(calc,1)',G_Ta(calc,4)',nextS)]
+
+ %overFlowTime = neville2(G_Ta(size(G_Ta,1)+(-1:0),1)',sum(G_Ta(size(G_Ta,1)+(-1:0),2:4),2)',4500)
+ end
+
+ A_stepAnIso(mu,type,eta,eps);
+ usedTime = G_Ta(end,:)
+ data = G_Da(end,:)
+
+ A_saveMesh(['anIso' int2str(type) '_' int2str(size(G_Ta,1))]);
+
+end
+
--- /dev/null
+function A_Iso(file,times,mu,type)
+%A_ANISO Summary of this function goes here
+% DeTiiled explanation goes here
+
+A_loadMeshF(file)
+
+global G_Ti;
+global G_Di;
+
+for i = 1:times
+ if(size(G_Ti,1)>2)
+ nextF = G_Ti(end,1)*4;
+ calc = size(G_Ti,1)+(-1:0);
+ nextTime = [nextS neville2(G_Ti(calc,1)',G_Ti(calc,2)',nextF) ...
+ neville2(G_Ti(calc,1)',G_Ti(calc,3)',nextF) ...
+ neville2(G_Ti(calc,1)',G_Ti(calc,4)',nextF)]
+
+ %overFlowTime = neville2(G_Ti(size(G_Ti,1)+(-1:0),1)',sum(G_Ti(size(G_Ti,1)+(-1:0),2:4),2)',4500)
+ end
+
+ A_stepIso(mu,type);
+ usedTime = G_Ti(end,:)
+ DiTi = G_Di(end,:)
+
+ A_saveMesh(['Iso' int2str(type) '_' int2str(size(G_Ti,1))]);
+
+end
+
-function A_loadMesh(coordinates,elements, neigh)
+function A_loadMesh(coordinates,elements,neigh)
% Laed ein Mesh
%
% A_loadMesh(coordinates,elements)
global G_C;
global G_E;
global G_N;
+ global G_Ta;
+ global G_Da;
+ global G_Ti;
+ global G_Di;
+
G_C = coordinates;
G_E = elements;
G_N = neigh;
+ G_Ta = [];
+ G_Da = [];
+ G_Ti = [];
+ G_Di = [];
end
-
--- /dev/null
+function A_loadMeshF(file)
+% Laed ein Mesh
+%
+% A_loadMesh(coordinates,elements)
+
+
+ load(file)
+
+ global G_C;
+ global G_E;
+ global G_N;
+ global G_Ta;
+ global G_Da;
+ global G_Ti;
+ global G_Di;
+ G_C = coordinates;
+ G_E = elements;
+ G_N = neigh;
+ G_Ta = [];
+ G_Da = [];
+ G_Ti = [];
+ G_Di = [];
+
+ if(exist('timeA','var')~=0)
+ G_Ta = timeA;
+ end
+ if(exist('dataA','var')~=0)
+ G_Da = dataA;
+ end
+ if(exist('timeI','var')~=0)
+ G_Ti = timeI;
+ end
+ if(exist('dataI','var')~=0)
+ G_Di = dataI;
+ end
+end
+
+
+
--- /dev/null
+function A_saveMesh(name)
+
+global G_E;
+global G_C;
+global G_N;
+global G_Ta;
+global G_Da;
+global G_Ti;
+global G_Di;
+
+timeI =[];
+dataI =[];
+timeA =[];
+dataA =[];
+coordinates = G_C;
+elements = G_E;
+neigh = G_N;
+if(~isempty(G_Ta))
+ timeA = G_Ta;
+ dataA = G_Da;
+end
+if(~isempty(G_Ti))
+ timeI = G_Ti;
+ dataI = G_Di;
+end
+
+save (['meshSave/' name], 'coordinates', 'elements','neigh','timeA','dataA','timeI','dataI')
+
+end
+
-function dataAniso = A_stepAnIso(mu,type)
+function [dataAniso time er] = A_stepAnIso(mu,type,eta,eps)
+% dataAniso = A_stepAnIso(mu,type,eta,eps)
+%
global G_E;
global G_C;
global G_N;
+global G_Ta;
+global G_Da;
if (isempty(G_E) || isempty(G_C) || isempty(G_N))
disp 'Error: Please use A_loadMesh first'
return
end
- [coordinates_fine,elements_fine,neigh_fine,f2s]=refineQuad(G_C,G_E,G_N,2*ones(1,size(G_E,1)));
+time = zeros(1,3);
+ tic
+ [coordinates_fine,elements_fine,neigh_fine,f2s]=refineQuad(G_C,G_E,G_N,2);
+ time(1) = toc;
+
+ tic
A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type);
+ time(2) = toc;
b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
- marked = mark(x_fine(f2s)',ind,0.5,0.5);
+ marked = mark(x_fine(f2s)',ind,eta,eps);
+
+% clear 'coordinates_fine' 'elements_fine' 'neigh_fine' 'f2s'
dataAniso = [size(G_E,1) sqrt(sum(ind)) xe_fine];
- [G_C, G_E, G_N] = refineQuad(G_C,G_E,G_N,marked);
-
+
+ if(size(G_E,1)~=length(marked))
+ disp 'Error: MarkierungsVektor ist fehlerhaft'
+ return
+ end
+
+ tic
+ [G_C, G_E, G_N, f, er] = refineQuad(G_C,G_E,G_N,marked);
+ time(3) = toc;
+% if(~isempty(er))
+% plotShape(G_C,G_E,'');
+% plotMark(er,G_C,G_E);
+% end
+
+
+ G_Ta(size(G_Ta,1)+1,1:4) = [size(G_E,1) time];
+ G_Da(size(G_Da,1)+1,1:3) = dataAniso;
+
+
+% plotShape(G_C,G_E,'');
end
-function dataIso = A_stepIso(mu,type)
+function [dataIso timer] = A_stepIso(mu,type)
global G_E;
global G_C;
global G_N;
+global G_Ti;
+global G_Di;
if (isempty(G_E) || isempty(G_C))
disp 'Error: Please use A_loadMesh first'
return
end
- [coordinates_fine,elements_fine,neigh_fine,f2s]=refineQuad(G_C,G_E,G_N,2);
+time = zeros(1,2);
+ tic
+ [coordinates_fine,elements_fine,neigh_fine,f2s,er]=refineQuad(G_C,G_E,G_N,2);
+ time(1) = toc;
+
+ tic
A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type);
+ time(2) = toc;
b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
dataIso = [size(G_E,1) sqrt(sum(ind)) xe_fine];
- [G_C, G_E, G_N] = refineQuad(G_C,G_E,G_N,2);
+
+ G_C = coordinates_fine; G_E = elements_fine; G_N = neigh_fine;
+
+% if(~isempty(er))
+% plotShape(G_C,G_E,'');
+% plotMark(er,G_C,G_E);
+% end
+
+
+ G_Ti(size(G_Ti,1)+1,1:3) = [size(G_E,1) time];
+ G_Di(size(G_Di,1)+1,1:3) = dataIso;
end
--- /dev/null
+function [ y ] = neville2(x,y,t )\r
+%[ p ] = neville2(x,y,t )\r
+\r
+\r
+ n=length(y);\r
+ t(2:n) = t(1);\r
+ \r
+ for m = 2:n\r
+ \r
+ j=n-m+1;\r
+ %... = zeigt MATLAB das die Zeile noch weiter geht\r
+ y = ((t(1:j) - x(1:j)) .* y(2:j+1) ... \r
+ -(t(1:j) - x(m:n)) .* y(1:j)) ...\r
+ ./ (x(m:n)-x(1:j));\r
+ \r
+ end\r
+\r
+\r
+end
\ No newline at end of file
--- /dev/null
+function plotMark(ele,coordinates,elements)
+
+hold on;
+ disp 'Plot Updated'
+for idx = ele
+ current = sum(coordinates(elements(idx,[2,4])',:),1)/2;
+ scatter3(current(1),current(2),current(3),'xg');
+end
+hold off;
+
+end
+
-function [coo,ele,nei,f2s] = refineQuad(coordinates,elements,neigh,type)
+function [coo,ele,nei,f2s,err] = refineQuad(coordinates,elements,neigh,type)
%
% [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type)
%
%
% P. Schaefer
+err = [];
%Type wenn nur ein Wert: aufblaehen
if([1 1] == size(type))
type = repmat(type, size(elements,1),1);
end
+if(size(elements,1)~=size(neigh,1)||size(elements,1)~=length(type))
+ disp 'Dimensionen passen nicht!'
+ coo = coordinates; ele = elements; nei = neigh; f2s = [];
+ return;
+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;
%Elementanzahl speichern
c_loop = size(elements,1);
G_ref_E = elements;
G_ref_C = coordinates;
G_ref_N = neigh;
+% G_ref_s = quadNorm(G_ref_C,G_ref_E,'w');
+% G_ref_s = sqrt(sum(G_ref_s.*G_ref_s,2))
+
+
G_ref_t = type;
G_ref_f2s = repmat([1:c_loop]',1,4);
%Parameter Freigeben (Speicher...)
clear elements coordinates neigh type
-%Jedes Element verfeinern
+% Flächeninhalt aufbauen
+updateS(1:c_loop);
-ea = 1:c_loop;
+ref_old = [];
-for i = ea(G_ref_t>1)
- refine(i);
+%Jedes Element verfeinern
+while(1==1)
+ ref = find(G_ref_t>1);
+ ref = reshape(ref,1,length(ref));
+ if(isequal(ref,ref_old))
+ disp([' Warning: ' int2str(length(ref)) ' Elements couldn t be refined'])
+ err = ref;
+ break;
+ end
+ ref_old = ref;
+ if(isempty(ref))
+ break;
+ end
+ for ele = ref(randperm(length(ref))) %LOL Zufall bringts
+ % # HangingNode Check
+ N = G_ref_N(ele,G_ref_N(ele,5:8)==0);
+ N2 = N(N~=0); %Nachbarn der Kanten mit nur einem Nachbar
+ if(~isempty(N2))
+ N3=mod(find((G_ref_N(N2',:)==ele)')-1,4)+5; % ACHTUNG noch mal überprüfen
+ N4=N2(diag(G_ref_N(N2',N3'))~=0);
+ if(~isempty(N4))
+ for i = N4
+ if(i>length(G_ref_t))
+ continue;
+ end
+
+ if(G_ref_t(i)<1)
+ continue;
+ end
+
+ G_ref_t(i) = 2;
+ end
+ continue;
+ end
+ end
+
+ % # Size Check
+% N5=N2(G_ref_s(N2')>G_ref_s(ele))
+% if(~isempty(N5))
+% %Setze N5 zum Verfeinern
+% for i = N5'
+% if(G_ref_t(i)<1)
+% continue;
+% else
+% G_ref_t(i) = 2;
+% end
+% end
+%
+% continue;
+% end
+
+
+ % Wenn Alle Überprüfungen durchgelaufen sind
+ refineE(ele);
+ updateN(ele);
+ end
end
-
%Rueckgabe zuweisen
coo = G_ref_C;
ele = G_ref_E;
%Doppelte Koordinaten loeschen
[coo l pos] = unique(coo,'rows');
+
ele = pos(ele);
%Globale Variablen freigeben
-clear G_ref_E G_ref_C G_ref_N G_ref_f2s G_ref_t
+clear G_ref_E G_ref_C G_ref_N G_ref_f2s G_ref_t G_ref_s
end
-
+%% Rekursive Funktion
function err = refine(ele)
err = 0;
% global G_ref_E;
return;
end
% if(G_ref_t(N4(i))<=2)
- G_ref_t(N4(i))=2; %WIRD GNADENLOS AUF 2 GESETZT
- test = refine(N4(i));
+% G_ref_t(N4(i))=2; %WIRD GNADENLOS AUF 2 GESETZT
+% test = refine(N4(i));
if(test==1)
G_ref_t(N4(i))=0;
err = 1;
% verfeinere dieses element
refineE(ele);
-
% setze Nachbarschafts relationen
updateN(ele);
end
+%% Element Verfeinern ! sollte nur ausgeführt werden wenn wirklich möglich
function refineE(ele)
% Element wird gnadenlos Verfeinert
global G_ref_E;
G_ref_f2s(ele,1:3)=c_ele+3:-1:c_ele+1;
% G_ref_f2s(ele,:)=[c_ele+3:-1:c_ele+1 ele];
+ updateS([ele,c_ele+1:c_ele+3]);
elseif(G_ref_t(ele)==3)
G_ref_C(c_coo+1,:) = (G_ref_C(el(1),:)+G_ref_C(el(4),:))/2;
G_ref_C(c_coo+2,:) = (G_ref_C(el(2),:)+G_ref_C(el(3),:))/2;
G_ref_E(c_ele+1,:) = [el(1),el(2),c_coo+2,c_coo+1];
G_ref_f2s(ele,[1 2])=c_ele+1;
% G_ref_f2s(ele,[3 4])=ele;
+ updateS([ele,c_ele+1]);
elseif(G_ref_t(ele)==4)
G_ref_C(c_coo+1,:) = (G_ref_C(el(1),:)+G_ref_C(el(2),:))/2;
G_ref_C(c_coo+2,:) = (G_ref_C(el(4),:)+G_ref_C(el(3),:))/2;
G_ref_E(c_ele+1,:) = [c_coo+1,el(2),el(3),c_coo+2];
G_ref_f2s(ele,[2 3])=c_ele+1;
% G_ref_f2s(ele,[1 4])=ele;
+ updateS([ele,c_ele+1]);
end
G_ref_t(ele) = 0;
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_C;
+%global G_ref_C;
global G_ref_N;
global G_ref_f2s;
end
-
-
-
end
+
+%% Berechnet die Flächeninhalte der Elemente
+function updateS(ele)
+%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
+ tri = G_ref_E(i,:);
+ a = (G_ref_C(tri(2),:)-G_ref_C(tri(1),:));
+ b = (G_ref_C(tri(4),:)-G_ref_C(tri(1),:));
+ N = cross(a',b');
+ N = N/norm(N);
+ G_ref_s(i) = sqrt(sum(N.*N));
+ end
+end
\ No newline at end of file
projects / bacc.git / commitdiff