From da86b7b5514e634b10cc83bfd66da8246194d891 Mon Sep 17 00:00:00 2001
From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Sat, 14 Jan 2012 22:40:47 +0000
Subject: [PATCH] =?utf8?q?[src]=20A=5F=20skripte=20reduziert=20[src]=20tes?=
 =?utf8?q?t=5Fsol=20hinzugef=C3=BCgt=20(10=20Schritte)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=utf8
Content-Transfer-Encoding: 8bit

git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@73 26120e32-c555-405d-b3e1-1f783fb42516
---
 src/A_AnIso.m                   | 36 -------------------
 src/A_Iso.m                     | 36 -------------------
 src/A_loadMesh.m                | 13 +++----
 src/A_loadMeshF.m               | 27 +++++---------
 src/A_plot.m                    | 45 +++++------------------
 src/A_run.m                     | 43 ++++++++++++++++++++++
 src/A_saveMesh.m                | 25 +++++--------
 src/{A_stepAnIso.m => A_step.m} | 23 +++++++-----
 src/A_stepIso.m                 | 63 ---------------------------------
 src/mark.m                      | 34 ++++++++++--------
 src/test_sol.m                  |  3 ++
 11 files changed, 111 insertions(+), 237 deletions(-)
 delete mode 100644 src/A_AnIso.m
 delete mode 100644 src/A_Iso.m
 create mode 100644 src/A_run.m
 rename src/{A_stepAnIso.m => A_step.m} (70%)
 delete mode 100644 src/A_stepIso.m
 create mode 100644 src/test_sol.m

diff --git a/src/A_AnIso.m b/src/A_AnIso.m
deleted file mode 100644
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
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 95f6fd0..1a34fc0 100644
--- 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 c85d24f..fa59a80 100644
--- 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 9789e21..72a1dc0 100644
--- 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
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 9652ced..4dce79e 100644
--- 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_stepAnIso.m b/src/A_step.m
similarity index 70%
rename from src/A_stepAnIso.m
rename to src/A_step.m
index 35d58aa..0b20dd4 100644
--- a/src/A_stepAnIso.m
+++ b/src/A_step.m
@@ -1,11 +1,11 @@
-function [dataAniso time er] = A_stepAnIso(mu,type,eta,eps)
+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_Ta;
-global G_Da;
+global G_T;
+global G_D;
 
 if (isempty(G_E) || isempty(G_C) || isempty(G_N))
   disp 'Error: Please use A_loadMesh first'
@@ -21,7 +21,6 @@ if(length(mu)==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;
@@ -29,9 +28,17 @@ time = zeros(1,3);
   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
 
   tic
-  dataAniso = size(G_E,1);
+  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;
 
@@ -39,7 +46,7 @@ time = zeros(1,3);
 
     ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
 
-    dataAniso = [dataAniso type(i) sqrt(sum(ind)) xe_fine];
+    data = [data type(i) sqrt(sum(ind)) xe_fine];
   end
   time(2) = toc;
   
@@ -64,8 +71,8 @@ time = zeros(1,3);
 %   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;
+  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,'');
diff --git a/src/A_stepIso.m b/src/A_stepIso.m
deleted file mode 100644
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 8c0c0d3..aae5d94 100644
--- 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
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
-- 
2.47.3