From b6528c4aba5181e4d23c360cdb1856cb7c9285a5 Mon Sep 17 00:00:00 2001 From: treecity Date: Sat, 22 Oct 2011 12:37:57 +0000 Subject: [PATCH] [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 --- src/A_loadMesh.m | 4 +- src/A_stepAnIso.m | 7 +- src/A_stepIso.m | 5 +- src/exmpl_2DLShape.mat | Bin 264 -> 321 bytes src/exmpl_2DQuad.mat | Bin 279 -> 331 bytes src/exmpl_2DQuad2.mat | Bin 0 -> 296 bytes src/exmpl_3DCube.mat | Bin 271 -> 339 bytes src/exmpl_3DFichCube.mat | Bin 322 -> 430 bytes src/nMesh/refineQuad.m | 1 + src/refineQuad.m | 309 ++++++++++++++++++++++++++++++++++----- 10 files changed, 281 insertions(+), 45 deletions(-) create mode 100644 src/exmpl_2DQuad2.mat diff --git a/src/A_loadMesh.m b/src/A_loadMesh.m index 703edb4..f1e48cf 100644 --- 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 5e4d976..cbd8182 100644 --- 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 5175ffb..dc0c917 100644 --- 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 GIT binary patch delta 98 zcmeBRI>($ diff --git a/src/exmpl_2DQuad.mat b/src/exmpl_2DQuad.mat index 743b64c55901799fa87238cadf110e05dfb327b8..4e1a628d29c3fa072730f65499808d946895eeb3 100644 GIT binary patch delta 159 zcmbQvbed^`i9~Q>iGqJ}iGq=lf}x?6p{139g@TcRq2a_p?THDj6KgmaT_?60ODV9( z-AG<;|6swzuYp>Ab_@)&Y=t-SGcYh%O}woZ%wQ<2&8@7>Zp|&N&aQ1Ntu4&Z5yRpE zk~ab3iaC#y6B0TYlFSsIHC$3S!1Z*qL#Eq`h`6k@Dz*t}Oi8R0Y>a2{DoIL8L`qy> KWSHT?bQ=HzSS`K) delta 107 zcmV-x0F?jB0+#}iG#FQ9WFSg)ZXh5yATc>QFfuwaI3O}GF)@)*E0JIYk$eaNPLYK( zBR~fq?6LD7*kk9z06HICKf1ZN%)@04=>PzfW$ai89CXQSlLSQtQZO{MGBmR?F;XxxFf>#k5il@%`tma{FxUWb#hk~<2?-nwNoESq z8ZIdu;Ci~5!O&R1n90ywK>hi1@uexxRHrQwYrbm8XvohT!5F})tQ=y^%y2G~)d6IN zF5HX_NM-=Fb5CLW7+f4w9Q@gLUjCt%~AD`!e3RH-SLF(jqAskx=8g^8JknTe$d96ZG$mUMuT J;iew54FJh9B1`}P delta 36 qcmcc2)Xy}*Lc}+*Qo+bj!O+wSh)m5Dj0_A7CkARyOkkZ@GZ_H3^a?-# diff --git a/src/exmpl_3DFichCube.mat b/src/exmpl_3DFichCube.mat index 3c5c502061df8d68352be66b8e2f9629e9783250..af813975bb31121987a76da10f696534bda5c8c1 100644 GIT binary patch delta 148 zcmV;F0Birk0{fNv_P?!4jMV<+zIeZ3yf zOiDj#LBh(?w2!&WYa49l?h(P#S*M9>kNeX}1s67=QA&+4_@lysv~l#l&L4iX69u-p C**_2f delta 39 tcmZ3-e28g+sd#Cgf)@}gm?{_85kN)4Ah>Oz&f$!7y#ag3wZzl diff --git a/src/nMesh/refineQuad.m b/src/nMesh/refineQuad.m index 4fde40a..c212b4f 100644 --- 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 f5f68cb..bf682f4 100644 --- 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 -- 2.47.3