[src] refineQuad REKURSIV fertig gestellt

[bug] mehr als 2000k ELemente gehen per Rekursion nicht
[src] Autoskripte an neue Nachbarversion angepasst

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

diff --git a/src/A_loadMesh.m b/src/A_loadMesh.m
index 703edb4417a60d1ff05a178a524a910a9e485077..f1e48cf4048fdb8cfdb99c1f8291b3f2d9fe3dfb 100644 (file)
--- a/src/A_loadMesh.m
+++ b/src/A_loadMesh.m
@@ -1,11 +1,13 @@
-function A_loadMesh(coordinates,elements)
+function A_loadMesh(coordinates,elements, neigh)
 % Laed ein Mesh
 %
 % A_loadMesh(coordinates,elements)
 
   global G_C;
   global G_E;
+  global G_N;
   G_C = coordinates;
   G_E = elements;
+  G_N = neigh;
 end
 

diff --git a/src/A_stepAnIso.m b/src/A_stepAnIso.m
index 5e4d97632aa7f14ef0d803ad193d4ce25b55e540..cbd818275c41b9dd2997f929ce5897c321f887c2 100644 (file)
--- a/src/A_stepAnIso.m
+++ b/src/A_stepAnIso.m
@@ -2,13 +2,14 @@ function dataAniso = A_stepAnIso(mu,type)
 
 global G_E;
 global G_C;
+global G_N;
 
-if (isempty(G_E) || isempty(G_C))
+if (isempty(G_E) || isempty(G_C) || isempty(G_N))
   disp 'Error: Please use A_loadMesh first'
   return
 end
 
-  [coordinates_fine,elements_fine, f2s]=refineQuad(G_C,G_E,2*ones(1,size(G_E,1)));
+  [coordinates_fine,elements_fine,neigh_fine,f2s]=refineQuad(G_C,G_E,G_N,2*ones(1,size(G_E,1)));
   
   A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type);
   
@@ -23,7 +24,7 @@ end
   marked = mark(x_fine(f2s)',ind,0.5,0.5);
   
   dataAniso = [size(G_E,1) sqrt(sum(ind)) xe_fine];
-  [G_C, G_E] = refineQuad(G_C,G_E,marked);
+  [G_C, G_E, G_N] = refineQuad(G_C,G_E,G_N,marked);
 
 
 end

diff --git a/src/A_stepIso.m b/src/A_stepIso.m
index 5175ffbafa69d508d1685262c689d9317b764bf6..dc0c917ea5c0c4c8fd04a16859bb230a223496f0 100644 (file)
--- a/src/A_stepIso.m
+++ b/src/A_stepIso.m
@@ -2,13 +2,14 @@ function dataIso = A_stepIso(mu,type)
 
 global G_E;
 global G_C;
+global G_N;
 
 if (isempty(G_E) || isempty(G_C))
   disp 'Error: Please use A_loadMesh first'
   return
 end
 
-  [coordinates_fine,elements_fine, f2s]=refineQuad(G_C,G_E,2);
+  [coordinates_fine,elements_fine,neigh_fine,f2s]=refineQuad(G_C,G_E,G_N,2);
   
   A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type);
   
@@ -21,7 +22,7 @@ end
   ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
   
   dataIso = [size(G_E,1) sqrt(sum(ind)) xe_fine];
-  [G_C, G_E] = refineQuad(G_C,G_E,2);
+  [G_C, G_E, G_N] = refineQuad(G_C,G_E,G_N,2);
 
 
 end

diff --git a/src/exmpl_2DLShape.mat b/src/exmpl_2DLShape.mat
index 9c9f86097cb660f6456cc5604866caa0ad92b34f..3b6e437c1684616b4d5778e43a77fda285659d1e 100644 (file)
Binary files a/src/exmpl_2DLShape.mat and b/src/exmpl_2DLShape.mat differ

diff --git a/src/exmpl_2DQuad.mat b/src/exmpl_2DQuad.mat
index 743b64c55901799fa87238cadf110e05dfb327b8..4e1a628d29c3fa072730f65499808d946895eeb3 100644 (file)
Binary files a/src/exmpl_2DQuad.mat and b/src/exmpl_2DQuad.mat differ

diff --git a/src/exmpl_2DQuad2.mat b/src/exmpl_2DQuad2.mat
new file mode 100644 (file)
index 0000000..eedc589
Binary files /dev/null and b/src/exmpl_2DQuad2.mat differ

diff --git a/src/exmpl_3DCube.mat b/src/exmpl_3DCube.mat
index d2e82344d2a377bd00f5c8a196cb38fe44ee98e0..448514ea32294854065929ee455be9d2ed664436 100644 (file)
Binary files a/src/exmpl_3DCube.mat and b/src/exmpl_3DCube.mat differ

diff --git a/src/exmpl_3DFichCube.mat b/src/exmpl_3DFichCube.mat
index 3c5c502061df8d68352be66b8e2f9629e9783250..af813975bb31121987a76da10f696534bda5c8c1 100644 (file)
Binary files a/src/exmpl_3DFichCube.mat and b/src/exmpl_3DFichCube.mat differ

diff --git a/src/nMesh/refineQuad.m b/src/nMesh/refineQuad.m
index 4fde40a74c730b1e9c6210b0f652bdc7212a3a21..c212b4f8ec24982f922273770ff89f3bfa3379d9 100644 (file)
--- a/src/nMesh/refineQuad.m
+++ b/src/nMesh/refineQuad.m
@@ -170,6 +170,7 @@ S = find(mod((this(1:4)~=0).*(this(5:8)==0),2))' %An welchen Kanten habe ich Nac
 D = find(this(5:8)~=0)'
 MS = mod(find((G_ref_N(this(S),:)==ele)')-1,4)+1  %An welchen Kanten bin ich Nachbar
 MD = mod(find((G_ref_N(this([D D+4]),:)==ele)')-1,4)+1
+MD = reshape(MD,2,lenght(MD)/2);
 
 G_ref_N(split,1:8) = 0;
 

diff --git a/src/refineQuad.m b/src/refineQuad.m
index f5f68cbae4bcde1937e8dd7fc467f87a71463afd..bf682f4458f19a81f530d3edba336ee56501d01b 100644 (file)
--- a/src/refineQuad.m
+++ b/src/refineQuad.m
@@ -1,10 +1,10 @@
-function [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type)
+function [coo,ele,nei,f2s] = refineQuad(coordinates,elements,neigh,type)
 %
 %  [coordinates,elements,fa2so] = refineQuad(coordinates,elements,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 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 sein.
 %
 % der Typ zu jedem Element entspricht dabei:
 % 1- keine Verfeinerung
@@ -14,59 +14,290 @@ function [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type)
 %
 % P. Schaefer
 
+
+%Type wenn nur ein Wert: aufblaehen
 if([1 1] == size(type))
     type = repmat(type, size(elements,1),1);
 end
 
+%Globale Variabelen aufbauen
+global G_ref_E;
+global G_ref_C;
+global G_ref_N;
+global G_ref_f2s;
+global G_ref_t;
+
+%Elementanzahl speichern
 c_loop = size(elements,1);
 
-fa2so = repmat([1:c_loop]',1,4);
+%Globale Variablen zuweisen
+G_ref_E = elements;
+G_ref_C = coordinates;
+G_ref_N = neigh;
+G_ref_t = type;
+G_ref_f2s = repmat([1:c_loop]',1,4);
+
+%Parameter Freigeben (Speicher...)
+clear elements coordinates neigh type
+
+%Jedes Element verfeinern
+
+ea = 1:c_loop;
+
+for i = ea(G_ref_t>1)
+  refine(i);
+end
+
+%Rueckgabe zuweisen
+coo = G_ref_C;
+ele = G_ref_E;
+nei = G_ref_N;
+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
+end
+
+
+function err = refine(ele)
+err = 0;
+% global G_ref_E;
+% global G_ref_C;
+global G_ref_N;
+% global G_ref_f2s;
+global G_ref_t;
+
+if(G_ref_t(ele)<2)
+  if(G_ref_t(ele)<1)
+    err = 1;
+  end
+  return;
+end
+
+% Ueberpruefe Nachbarn auf Nodes
+N = G_ref_N(ele,G_ref_N(ele,5:8)==0);
+N2 = N(N~=0);
+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);
+  
+  % wenn ungueltig verfeinere sie (link auf Soll verfeinert werden?)
+  if(~isempty(N4))
+    for i = 1:length(N4)
+        if(G_ref_t(N4(i))==0)
+          err = 1;
+          return;
+        end
+%       if(G_ref_t(N4(i))<=2)
+        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;
+          return;
+        end
+%       else
+%         mod(find(G_ref_N(N4(i),:)==ele)-1,2)
+%       end
+    end
+  end
+end
+
+if((G_ref_t(ele)<1) || err==1)
+  err = 1;
+  return;
+end
+
+% verfeinere dieses element
+refineE(ele);
+
+% setze Nachbarschafts relationen
+updateN(ele);
+
+end
+
+function refineE(ele)
+% Element wird gnadenlos Verfeinert
+global G_ref_E;
+global G_ref_C;
+%global G_ref_N;
+global G_ref_f2s;
+global G_ref_t;
 
-for i = 1:c_loop
-    c_ele = size(elements,1);
-    c_coo = size(coordinates,1);
+    c_ele = size(G_ref_E,1);
+    c_coo = size(G_ref_C,1);
     
     
-    el = elements(i,:);
-    if(type(i)==1)
-        continue;
-    elseif(type(i)==2)
-        coordinates(c_coo+1,:) = (coordinates(el(1),:)+coordinates(el(2),:))/2;
-        coordinates(c_coo+2,:) = (coordinates(el(2),:)+coordinates(el(3),:))/2;
-        coordinates(c_coo+3,:) = (coordinates(el(3),:)+coordinates(el(4),:))/2;
-        coordinates(c_coo+4,:) = (coordinates(el(1),:)+coordinates(el(4),:))/2;
-        coordinates(c_coo+5,:) = (coordinates(el(1),:)+coordinates(el(3),:))/2;
-
-        elements(i,:) = [c_coo+4,c_coo+5,c_coo+3,el(4)];
-        elements(c_ele+1,:) = [c_coo+5,c_coo+2,el(3),c_coo+3];
-        elements(c_ele+2,:) = [c_coo+1,el(2),c_coo+2,c_coo+5];
-        elements(c_ele+3,:) = [el(1),c_coo+1,c_coo+5,c_coo+4];
+    el = G_ref_E(ele,:);
+    if(G_ref_t(ele)==1)
+      disp 'Warning: Type = 1';
+        return;
+    elseif(G_ref_t(ele)==2)
+        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(2),:)+G_ref_C(el(3),:))/2;
+        G_ref_C(c_coo+3,:) = (G_ref_C(el(3),:)+G_ref_C(el(4),:))/2;
+        G_ref_C(c_coo+4,:) = (G_ref_C(el(1),:)+G_ref_C(el(4),:))/2;
+        G_ref_C(c_coo+5,:) = (G_ref_C(el(1),:)+G_ref_C(el(3),:))/2;
+
+        G_ref_E(ele,:) = [c_coo+4,c_coo+5,c_coo+3,el(4)];
+        G_ref_E(c_ele+1,:) = [c_coo+5,c_coo+2,el(3),c_coo+3];
+        G_ref_E(c_ele+2,:) = [c_coo+1,el(2),c_coo+2,c_coo+5];
+        G_ref_E(c_ele+3,:) = [el(1),c_coo+1,c_coo+5,c_coo+4];
     
-        fa2so(i,2:4)=c_ele+1:c_ele+3;
-    elseif(type(i)==3)
-        coordinates(c_coo+1,:) = (coordinates(el(1),:)+coordinates(el(4),:))/2;
-        coordinates(c_coo+2,:) = (coordinates(el(2),:)+coordinates(el(3),:))/2;
+        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];
+    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;
         
-        elements(i,1) = c_coo+1;
-        elements(i,2) = c_coo+2;
+        G_ref_E(ele,1) = c_coo+1;
+        G_ref_E(ele,2) = c_coo+2;
         
-        elements(c_ele+1,:) = [el(1),el(2),c_coo+2,c_coo+1];
-        fa2so(i,[3 4])=c_ele+1;
-    elseif(type(i)==4)
-        coordinates(c_coo+1,:) = (coordinates(el(1),:)+coordinates(el(2),:))/2;
-        coordinates(c_coo+2,:) = (coordinates(el(4),:)+coordinates(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;
+     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;
         
-        elements(i,2) = c_coo+1;
-        elements(i,3) = c_coo+2;
+        G_ref_E(ele,2) = c_coo+1;
+        G_ref_E(ele,3) = c_coo+2;
         
-        elements(c_ele+1,:) = [c_coo+1,el(2),el(3),c_coo+2];
-        fa2so(i,[2 3])=c_ele+1;
+        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;
     end
+    
+    G_ref_t(ele) = 0;
+    
 end
 
-%igitigit
- [G_ref_C l pos] = unique(G_ref_C,'rows');
- G_ref_E = pos(G_ref_E);
- 
+function updateN(ele)
+% Nachbarschaften werden neu gesetzt (nach N und f2s)
+global G_ref_E;
+global G_ref_C;
+global G_ref_N;
+global G_ref_f2s;
+
+
+%Innere Beziehungen setzen
+% ...
+
+% plotShape(G_ref_C,G_ref_E);
+
+this = G_ref_N(ele,:);
+% ele
+split = G_ref_f2s(ele,:);
+S = find(mod((this(1:4)~=0).*(this(5:8)==0),2))'; %An welchen Kanten habe ich Nachbarn
+D = find(this(5:8)~=0)';
+MS = mod(find((G_ref_N(this(S),:)==ele)')-1,4)+1;  %An welchen Kanten bin ich Nachbar
+MD = mod(find((G_ref_N(this([D D+4]),:)==ele)')-1,4)+1;
+MD = reshape(MD,length(MD)/2,2);
+
+G_ref_N(split,1:8) = 0;
+
+if(split(1)==split(2))
+  G_ref_N(split([1 3])',1:4) = [ 0 0 split(3) 0;split(1) 0 0 0];
+      % Beziehungen fuer Kanten mit einem Nachbar
+    for i = 1:length(S)
+      if(mod(S(i),2)==0)
+        G_ref_N(this(S(i)),[MS(i) MS(i)+4]) = [split(S(i)) split(mod(S(i),4)+1)];
+        G_ref_N([split(S(i)) split(mod(S(i),4)+1)]',S(i))=this(S(i));
+      else
+        G_ref_N(this(S(i)),MS(i)) = split(S(i));
+        G_ref_N(split(S(i)),S(i)) = this(S(i));
+      end
+    end
+  
+    % Beziehungen fuer Kanten mit zwei Nachbarn
+    for i = 1:length(D)
+      if(mod(D(i),2)==0)
+        if(length(unique([G_ref_E(this(D(i)),:) G_ref_E(split(D(i)),:)]))==7)
+          G_ref_N(this(D(i)),MD(i,1)) = split(D(i));
+          G_ref_N(this(D(i)+4),MD(i,2)) = split(mod(D(i),4)+1);
+          G_ref_N(split(D(i)),D(i)) = this(D(i));
+          G_ref_N(split(mod(D(i),4)+1),D(i)) = this(D(i)+4);
+        else
+          G_ref_N(this(D(i)),MD(i,1)) = split(mod(D(i),4)+1);
+          G_ref_N(this(D(i)+4),MD(i,2)) = split(D(i));
+          G_ref_N(split(D(i)),D(i)) = this(D(i)+4);
+          G_ref_N(split(mod(D(i),4)+1),D(i)) = this(D(i));      
+        end
+      else
+        G_ref_N(this(D(i)),MD(i,1)) = split(D(i));
+        G_ref_N(this(D(i)+4),MD(i,2)) = split(D(i));
+        G_ref_N(split(D(i)),[D(i) D(i)+4]) = [this(D(i)) this(D(i)+4)];
+      end
+    end
+  
+elseif(split(1)==split(4))
+  G_ref_N(split([1 2])',1:4) = [ 0 split(2) 0 0;0 0 0 split(1)];
+    % Beziehungen fuer Kanten mit einem Nachbar
+    for i = 1:length(S)
+      if(mod(S(i),2)==1)
+        G_ref_N(this(S(i)),[MS(i) MS(i)+4]) = [split(S(i)) split(mod(S(i),4)+1)];
+        G_ref_N([split(S(i)) split(mod(S(i),4)+1)]',S(i))=this(S(i));
+      else
+        G_ref_N(this(S(i)),MS(i)) = split(S(i));
+        G_ref_N(split(S(i)),S(i)) = this(S(i));
+      end
+    end
+  
+    % Beziehungen fuer Kanten mit zwei Nachbarn
+    for i = 1:length(D)
+      if(mod(D(i),2)==1)
+        if(length(unique([G_ref_E(this(D(i)),:) G_ref_E(split(D(i)),:)]))==7)
+          G_ref_N(this(D(i)),MD(i,1)) = split(D(i));
+          G_ref_N(this(D(i)+4),MD(i,2)) = split(mod(D(i),4)+1);
+          G_ref_N(split(D(i)),D(i)) = this(D(i));
+          G_ref_N(split(mod(D(i),4)+1),D(i)) = this(D(i)+4);
+        else
+          G_ref_N(this(D(i)),MD(i,1)) = split(mod(D(i),4)+1);
+          G_ref_N(this(D(i)+4),MD(i,2)) = split(D(i));
+          G_ref_N(split(D(i)),D(i)) = this(D(i)+4);
+          G_ref_N(split(mod(D(i),4)+1),D(i)) = this(D(i));      
+        end
+      else
+        G_ref_N(this(D(i)),MD(i,1)) = split(D(i));
+        G_ref_N(this(D(i)+4),MD(i,2)) = split(D(i));
+        G_ref_N(split(D(i)),[D(i) D(i)+4]) = [this(D(i)) this(D(i)+4)];
+      end
+    end
+  
+else
+  G_ref_N(split',1:4) = [0 split(2) split(4) 0; 0 0 split(3) split(1);...
+                         split(2) 0 0 split(4); split(1) split(3) 0 0];
+  
+    % Beziehungen fuer Kanten mit einem Nachbar
+    for i = 1:length(S)
+      G_ref_N(this(S(i)),[MS(i) MS(i)+4]) = [split(S(i)) split(mod(S(i),4)+1)];
+      G_ref_N([split(S(i)) split(mod(S(i),4)+1)]',S(i))=this(S(i));
+    end
+  
+    % Beziehungen fuer Kanten mit zwei Nachbarn
+    for i = 1:length(D)
+      if(length(unique([G_ref_E(this(D(i)),:) G_ref_E(split(D(i)),:)]))==7)
+        G_ref_N(this(D(i)),MD(i,1)) = split(D(i));
+        G_ref_N(this(D(i)+4),MD(i,2)) = split(mod(D(i),4)+1);
+        G_ref_N(split(D(i)),D(i)) = this(D(i));
+        G_ref_N(split(mod(D(i),4)+1),D(i)) = this(D(i)+4);
+      else
+        G_ref_N(this(D(i)),MD(i,1)) = split(mod(D(i),4)+1);
+        G_ref_N(this(D(i)+4),MD(i,2)) = split(D(i));
+        G_ref_N(split(D(i)),D(i)) = this(D(i)+4);
+        G_ref_N(split(mod(D(i),4)+1),D(i)) = this(D(i));      
+      end
+    end
+  
+end
+
+
+
+
+
 end