[src] A_ skripte reduziert

[src] test_sol hinzugefügt (10 Schritte)

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

diff --git a/src/A_AnIso.m b/src/A_AnIso.m
deleted file mode 100644 (file)
index 74345cd..0000000
--- a/src/A_AnIso.m
+++ /dev/null
@@ -1,36 +0,0 @@
-function A_AnIso(file,times,mu,type,eta,eps)
-% A_AnIso(file,times,mu,type,eta,eps)
-
-
-A_loadMeshF(file)
-
-global G_Ta;
-global G_Da;
-
-if(length(mu)~=max(type)-1&&length(mu)~=1)
-  disp 'Error: Pleas set right type and mu parameters'
-  return
-end
-
-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);
-  cald = size(G_Ta,1)+(-1:0);
-  calc = size(G_Ta,1)+(-2:0);
-  nextTime = [nextS neville2(G_Ta(calc,1)',G_Ta(calc,2)',nextF) ...
-    neville2(G_Ta(cald,1)',G_Ta(cald,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,:)
- typeN = int2str(type);
-  
- A_saveMesh(['anIso' typeN(typeN~=' ') '_' int2str(size(G_Ta,1))]);
- 
-end
-

diff --git a/src/A_Iso.m b/src/A_Iso.m
deleted file mode 100644 (file)
index 8f78739..0000000
--- a/src/A_Iso.m
+++ /dev/null
@@ -1,36 +0,0 @@
-function A_Iso(file,times,mu,type)
-%A_Iso(file,times,mu,type)
-
-format longG
-
-A_loadMeshF(file)
-
-global G_Ti;
-global G_Di;
-
-if(length(mu)~=max(type)-1&&length(mu)~=1)
-  disp 'Error: Pleas set right type and mu parameters'
-  return
-end
-
-for i = 1:times
- if(size(G_Ti,1)>2)
-  nextF = G_Ti(end,1)*4;
-  cald = size(G_Ti,1)+(-1:0);
-  calc = size(G_Ti,1)+(-2:0);
-  nextTime = [nextF neville2(G_Ti(cald,1)',G_Ti(cald,2)',nextF) ...
-    neville2(G_Ti(calc,1)',G_Ti(calc,3)',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,:)
- 
- typeN = int2str(type);
-  
- A_saveMesh(['Iso' typeN(typeN~=' ') '_' int2str(size(G_Ti,1))]);
- 
-end
-

diff --git a/src/A_loadMesh.m b/src/A_loadMesh.m
index 95f6fd078dac75ca32257a949d64f4b287adf45f..1a34fc02bc1ee5bfd21e06b6e715ce249303ae71 100644 (file)
--- a/src/A_loadMesh.m
+++ b/src/A_loadMesh.m
@@ -6,16 +6,13 @@ 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;
+  global G_T;
+  global G_D;
   
   G_C = coordinates;
   G_E = elements;
   G_N = neigh;
-  G_Ta = [];
-  G_Da = [];
-  G_Ti = [];
-  G_Di = [];
+  G_T = [];
+  G_D = [];
+
 end

diff --git a/src/A_loadMeshF.m b/src/A_loadMeshF.m
index c85d24fe48a285b82398b6e8ea287ec1e0e2130b..fa59a8042c2868b665e052c26545cdc80f57fa2d 100644 (file)
--- a/src/A_loadMeshF.m
+++ b/src/A_loadMeshF.m
@@ -9,29 +9,20 @@ function A_loadMeshF(file)
   global G_C;
   global G_E;
   global G_N;
-  global G_Ta;
-  global G_Da;
-  global G_Ti;
-  global G_Di;
+  global G_T;
+  global G_D;
+  
   G_C = coordinates;
   G_E = elements;
   G_N = neigh;
-  G_Ta = [];
-  G_Da = []; 
-  G_Ti = [];
-  G_Di = [];
+  G_T = [];
+  G_D = [];
   
-  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;
+  if(exist('time','var')~=0)
+    G_T = time;
   end
-  if(exist('dataI','var')~=0)
-    G_Di = dataI;
+  if(exist('datA','var')~=0)
+    G_D = data;
   end
 end
 

diff --git a/src/A_plot.m b/src/A_plot.m
index 9789e21c0701fefeee10dd08a9597c0781dd7d67..72a1dc015fa15d352d1b79193eb4e782036448ec 100644 (file)
--- a/src/A_plot.m
+++ b/src/A_plot.m
@@ -2,39 +2,10 @@ function A_plot()
 
 type2str = ['Analytisch  ' ; 'Quad Element' ; 'Quad Achse  '; 'Quad Seite  '];
 
+%%G_D
+global G_D;
 
-%% G_Di
-global G_Di;
-
-[m n] = size(G_Di);
-
-step = round(n/3);
-
-
-if step<1
-    disp ('Error: No Data to show.')
-else
-
-figure(2)
-loglog(G_Di(:,1),G_Di(:,[3+(0:step-1)*3]))
-
-xlabel('Elemente');
-ylabel('Schaetzer');
-legend(type2str (G_Di(1,[2+(0:step-1)*3])',:));
-
-
-figure(3)
-loglog(G_Di(:,1),G_Di(:,[4+(0:step-1)*3]))
-
-xlabel('Elemente');
-ylabel('eNorm');
-legend(type2str(G_Di(1,[2+(0:step-1)*3])',:));
-end
-
-%%G_Da
-global G_Da;
-
-[m n] = size(G_Da);
+[m n] = size(G_D);
 
 step = round(n/3);
 
@@ -44,19 +15,19 @@ if step<1
 else
 
 figure(4)
-loglog(G_Da(:,1),G_Da(:,[3+(0:step-1)*3]))
+loglog(G_D(:,1),G_D(:,[3+(0:step-1)*3]))
 
 xlabel('Elemente');
 ylabel('Schaetzer');
-legend(type2str (G_Da(1,[2+(0:step-1)*3])',:));
+legend(type2str (G_D(1,[2+(0:step-1)*3])',:));
 
 
 figure(5)
-loglog(G_Da(:,1),G_Da(:,[4+(0:step-1)*3]))
+loglog(G_D(:,1),G_D(:,[4+(0:step-1)*3]))
 
-xlabel('Elemente');3
+xlabel('Elemente');
 ylabel('eNorm');
-legend(type2str(G_Da(1,[2+(0:step-1)*3])',:));
+legend(type2str(G_D(1,[2+(0:step-1)*3])',:));
 end
 
 end
\ No newline at end of file

diff --git a/src/A_run.m b/src/A_run.m
new file mode 100644 (file)
index 0000000..316697f
--- /dev/null
+++ b/src/A_run.m
@@ -0,0 +1,43 @@
+function A_run(file,times,mu,type,eta,eps,out)
+% A_AnIso(file,times,mu,type,eta,eps,out)
+% file - starting mesh
+% times - how often
+% mu - Zulaessigkeitsbed
+% type - art des Tests
+% eta - alle verfeinern oder nur wichtige? 
+% eps - wie verfeinern iso or aniso
+% out - Dateiname der Ausgabe
+
+
+A_loadMeshF(file)
+
+global G_T;
+global G_D;
+
+if(length(mu)~=max(type)-1&&length(mu)~=1)
+  disp 'Error: Pleas set right type and mu parameters'
+  return
+end
+
+for i = 1:times
+ 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);
+  cald = size(G_T,1)+(-1:0);
+  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)
+ end
+
+ A_step(mu,type,eta,eps);
+ usedTime = G_T(end,:)
+ data = G_D(end,:)
+ typeN = int2str(type);
+  
+ A_saveMesh([out typeN(typeN~=' ') '_' int2str(size(G_T,1))]);
+ 
+end
+

diff --git a/src/A_saveMesh.m b/src/A_saveMesh.m
index 9652cedc935316d4dafbb84b37ec66135ec11f63..4dce79eac34271eb6394610cb2671e4d4ad5df71 100644 (file)
--- a/src/A_saveMesh.m
+++ b/src/A_saveMesh.m
@@ -3,28 +3,21 @@ 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;
+global G_T;
+global G_D;
 
-timeI =[];
-dataI =[];
-timeA =[];
-dataA =[];
+time =[];
+data =[];
 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;
+if(~isempty(G_T))
+  time = G_T;
+  data = G_D;
 end
 
-save (['meshSave/' name], 'coordinates', 'elements','neigh','timeA','dataA','timeI','dataI')
+
+save (['meshSave/' name], 'coordinates', 'elements','neigh','time','data')
 
 end
 

diff --git a/src/A_step.m b/src/A_step.m
new file mode 100644 (file)
index 0000000..0b20dd4
--- /dev/null
+++ b/src/A_step.m
@@ -0,0 +1,81 @@
+function [data time er] = A_step(mu,type,eta,eps)
+% [dataAniso time er] = A_stepAnIso(mu,type,eta,eps)
+
+global G_E;
+global G_C;
+global G_N;
+global G_T;
+global G_D;
+
+if (isempty(G_E) || isempty(G_C) || isempty(G_N))
+  disp 'Error: Please use A_loadMesh first'
+  return
+end
+
+if(length(mu)~=max(type)-1&&length(mu)~=1)
+  disp 'Error: Pleas set right type and mu parameters'
+  return
+end
+if(length(mu)==1)
+  mu = repmat(mu,max(type)-1,1);
+end
+
+time = zeros(1,3);
+  tic
+  [coordinates_fine,elements_fine,neigh_fine,f2s]=refineQuad(G_C,G_E,G_N,2);
+  time(1) = toc;
+
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+
+  tic
+  data = size(G_E,1);
+  for i = 1:length(type)
+    A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type(i));
+    
+      [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
+
+    x_fine = A_fine\b;
+
+    xe_fine = x_fine'*A_fine*x_fine;
+
+    ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
+
+    data = [data type(i) sqrt(sum(ind)) xe_fine];
+  end
+  time(2) = toc;
+  
+  
+  marked = mark(x_fine(f2s)',ind,eta,eps);
+  
+%   clear 'coordinates_fine' 'elements_fine' 'neigh_fine' 'f2s'
+  
+
+  
+  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_T(size(G_T,1)+1,1:4) = [size(G_E,1) time];
+  G_D(size(G_D,1)+1,1:length(data)) = data;
+  
+  
+%   plotShape(G_C,G_E,'');
+
+end
+

diff --git a/src/A_stepAnIso.m b/src/A_stepAnIso.m
deleted file mode 100644 (file)
index 35d58aa..0000000
--- a/src/A_stepAnIso.m
+++ /dev/null
@@ -1,74 +0,0 @@
-function [dataAniso time er] = A_stepAnIso(mu,type,eta,eps)
-% [dataAniso time er] = 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
-
-if(length(mu)~=max(type)-1&&length(mu)~=1)
-  disp 'Error: Pleas set right type and mu parameters'
-  return
-end
-if(length(mu)==1)
-  mu = repmat(mu,max(type)-1,1);
-end
-
-time = zeros(1,3);
-  
-  tic
-  [coordinates_fine,elements_fine,neigh_fine,f2s]=refineQuad(G_C,G_E,G_N,2);
-  time(1) = toc;
-
-  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
-
-  tic
-  dataAniso = size(G_E,1);
-  for i = 1:length(type)
-    A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type(i));
-
-    x_fine = A_fine\b;
-
-    xe_fine = x_fine'*A_fine*x_fine;
-
-    ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
-
-    dataAniso = [dataAniso type(i) sqrt(sum(ind)) xe_fine];
-  end
-  time(2) = toc;
-  
-  
-  marked = mark(x_fine(f2s)',ind,eta,eps);
-  
-%   clear 'coordinates_fine' 'elements_fine' 'neigh_fine' 'f2s'
-  
-
-  
-  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:length(dataAniso)) = dataAniso;
-  
-  
-%   plotShape(G_C,G_E,'');
-
-end
-

diff --git a/src/A_stepIso.m b/src/A_stepIso.m
deleted file mode 100644 (file)
index 1abcfad..0000000
--- a/src/A_stepIso.m
+++ /dev/null
@@ -1,63 +0,0 @@
-function [dataIso timer] = A_stepIso(mu,type)
-% [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
-
-if(length(mu)~=max(type)-1&&length(mu)~=1)
-  disp 'Error: Pleas set right type and mu parameters'
-  return
-end
-if(length(mu)==1)
-  mu = repmat(mu,max(type)-1,1);
-end
-
-time = zeros(1,2);
-  tic
-  [coordinates_fine,elements_fine,neigh_fine,f2s,er]=refineQuad(G_C,G_E,G_N,2);
-  time(1) = toc;
-  
-  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2)); 
-  
-  tic
-  dataIso = size(G_E,1);
-  for i = 1:length(type)
-  A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type(i));
-  [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)
-  end
-  x_fine = A_fine\b;
-  
-  xe_fine = x_fine'*A_fine*x_fine;
-
-  ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
-  
-  dataIso = [dataIso type(i) sqrt(sum(ind)) xe_fine];
-  end
-  time(2) = toc;
-  
-  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:length(dataIso)) = dataIso;
-
-
-end
-

diff --git a/src/mark.m b/src/mark.m
index 8c0c0d393d3780b392fcb61b3c071f877709aca8..aae5d94e3079a128609750625848a9c15a0412ea 100644 (file)
--- a/src/mark.m
+++ b/src/mark.m
@@ -1,6 +1,9 @@
 function REF = mark(xF2S,ind,eta,eps)
-%REFINETYPE Summary of this function goes here
-%   Detailed explanation goes here
+% function REF = mark(xF2S,ind,eta,eps)
+% xF2S - Father son 
+% ind - error estimator
+% eta - refine element? (0..1, 1 = All)
+% eps - refine how? (0 = Isotrop)
 
 if(size(xF2S,1)==1)
   xF2S = xF2S';
@@ -10,29 +13,30 @@ T4 = [1 1 1 1; 1 -1 1 -1; 1 1 -1 -1;1 -1 -1 1]/4;
 
 REF=ones(1,size(xF2S,2));
 t1=zeros(1,size(xF2S,2));
+t3 =0; t4 = 0;
 
 Ct = T4*xF2S;
 
 %% Muss ueberhaupt verfeinert werden
+if(eta <1)
+    [s_ind idx] = sort(ind,'descend');
 
-[s_ind idx] = sort(ind,'descend');
+    sum_ind = cumsum(s_ind,1);
 
-sum_ind = cumsum(s_ind,1);
+    ell = find(sum_ind  >= sum_ind(end) * eta,1);
 
-ell = max(find(sum_ind >= sum_ind(end)*eta,1));
-
-t1(idx(ell:end)) = 1;
-
-% REF = eps * ind >
-
-% t1(7) = 1;
+    t1(idx(ell:end)) = 1;
+end
 
 
 %% Wie muss verfeinert werden
-t3 = (eps*abs(Ct(3,:)) >= sqrt(Ct(2,:).^2+Ct(4,:).^2));
-REF(t3) = 3;
-t4 = (eps*abs(Ct(4,:)) >= sqrt(Ct(2,:).^2+Ct(3,:).^2));
-REF(t4) = 4;
+if(eps > 0)
+    t3 = (eps*abs(Ct(3,:)) >= sqrt(Ct(2,:).^2+Ct(4,:).^2));
+    REF(t3) = 3;
+    t4 = (eps*abs(Ct(4,:)) >= sqrt(Ct(2,:).^2+Ct(3,:).^2));
+    REF(t4) = 4;
+end
 REF(~(t4+t3+t1)) = 2;
 
+
 end

diff --git a/src/test_sol.m b/src/test_sol.m
new file mode 100644 (file)
index 0000000..57ba42e
--- /dev/null
+++ b/src/test_sol.m
@@ -0,0 +1,3 @@
+mex mex_build_AU.cpp slpRectangle.cpp
+A_run('exmpl_2DQuad.mat',10,[0.7],[1 2 3],0.9,0,'AdaptIso_');
+A_run('exmpl_2DQuad.mat',10,[0.7],[1 2 3],0.9,0.7,'AdaptAnIso_');
\ No newline at end of file