From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Sat, 22 Oct 2011 21:12:47 +0000 (+0000)
Subject: [src] refineQuad komplett auf Iterativ umgestellt
X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=05cc84a4f9c822abcf090d41aecdd8e32c15e27f;p=bacc.git

[src] refineQuad komplett auf Iterativ umgestellt
[src] A_ skripte Hinzugefügt die Eine Datei laden und mit bestimmten Parameter mehrfach ausführen
[src] A_ Daten werden automatisch zwischen gespeichert -> können geladen werden und fortgesetzt
[src] plotMark -> Ergänzung zu plotShape: Markiert einzelne Elemente im Gitter
[src] neville zum kalkulieren hinzugefügt.. (eigentlich sinnlos weils einfacher gehen könnte...)

git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@59 26120e32-c555-405d-b3e1-1f783fb42516
---

diff --git a/src/A_AnIso.m b/src/A_AnIso.m
new file mode 100644
index 0000000..4fc8049
--- /dev/null
+++ b/src/A_AnIso.m
@@ -0,0 +1,29 @@
+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
+
diff --git a/src/A_Iso.m b/src/A_Iso.m
new file mode 100644
index 0000000..2385163
--- /dev/null
+++ b/src/A_Iso.m
@@ -0,0 +1,28 @@
+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
+
diff --git a/src/A_loadMesh.m b/src/A_loadMesh.m
index f1e48cf..95f6fd0 100644
--- a/src/A_loadMesh.m
+++ b/src/A_loadMesh.m
@@ -1,4 +1,4 @@
-function A_loadMesh(coordinates,elements, neigh)
+function A_loadMesh(coordinates,elements,neigh)
 % Laed ein Mesh
 %
 % A_loadMesh(coordinates,elements)
@@ -6,8 +6,16 @@ function A_loadMesh(coordinates,elements, neigh)
   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
-
diff --git a/src/A_loadMeshF.m b/src/A_loadMeshF.m
new file mode 100644
index 0000000..c85d24f
--- /dev/null
+++ b/src/A_loadMeshF.m
@@ -0,0 +1,39 @@
+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
+
+
+
diff --git a/src/A_saveMesh.m b/src/A_saveMesh.m
new file mode 100644
index 0000000..9652ced
--- /dev/null
+++ b/src/A_saveMesh.m
@@ -0,0 +1,30 @@
+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
+
diff --git a/src/A_stepAnIso.m b/src/A_stepAnIso.m
index cbd8182..8286ec4 100644
--- a/src/A_stepAnIso.m
+++ b/src/A_stepAnIso.m
@@ -1,17 +1,27 @@
-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));
   
@@ -21,11 +31,31 @@ end
 
   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
 
diff --git a/src/A_stepIso.m b/src/A_stepIso.m
index dc0c917..0eea0b8 100644
--- a/src/A_stepIso.m
+++ b/src/A_stepIso.m
@@ -1,17 +1,25 @@
-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));
   
@@ -22,7 +30,17 @@ end
   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
diff --git a/src/neville2.m b/src/neville2.m
new file mode 100644
index 0000000..b15b0e6
--- /dev/null
+++ b/src/neville2.m
@@ -0,0 +1,19 @@
+function [ y ] = neville2(x,y,t )
+%[ p ] = neville2(x,y,t )
+
+
+    n=length(y);
+    t(2:n) = t(1);
+    
+    for m = 2:n
+        
+        j=n-m+1;
+            %... = zeigt MATLAB das die Zeile noch weiter geht
+        y = ((t(1:j) - x(1:j)) .* y(2:j+1) ... 
+                -(t(1:j) - x(m:n)) .* y(1:j)) ...
+                ./ (x(m:n)-x(1:j));
+            
+    end
+
+
+end
\ No newline at end of file
diff --git a/src/plotMark.m b/src/plotMark.m
new file mode 100644
index 0000000..26eb9db
--- /dev/null
+++ b/src/plotMark.m
@@ -0,0 +1,12 @@
+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
+
diff --git a/src/refineQuad.m b/src/refineQuad.m
index bf682f4..9a2f7d9 100644
--- a/src/refineQuad.m
+++ b/src/refineQuad.m
@@ -1,4 +1,4 @@
-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)
 %
@@ -14,18 +14,26 @@ function [coo,ele,nei,f2s] = refineQuad(coordinates,elements,neigh,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);
@@ -34,20 +42,78 @@ 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;
@@ -56,13 +122,14 @@ f2s = G_ref_f2s;
 
 %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;
@@ -93,8 +160,8 @@ if(~isempty(N2))
           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;
@@ -114,12 +181,12 @@ end
 
 % 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;
@@ -150,6 +217,7 @@ global G_ref_t;
     
         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;
@@ -160,6 +228,7 @@ global G_ref_t;
         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;
@@ -170,16 +239,19 @@ global G_ref_t;
         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;
 
@@ -296,8 +368,25 @@ else
 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