From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Sun, 28 Aug 2011 15:28:57 +0000 (+0000)
Subject: [src] eigenartig -> Funktionen anders füttern geht aber Variablen tauschen nicht... 
X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=8345056d2ffc66adc5a730a69b9653c80c53a6b2;p=bacc.git

[src] eigenartig -> Funktionen anders füttern geht aber Variablen tauschen nicht. . .

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

diff --git a/src/mex_build_As1.cpp b/src/mex_build_As1.cpp
index d54e512..485b74f 100644
--- a/src/mex_build_As1.cpp
+++ b/src/mex_build_As1.cpp
@@ -69,7 +69,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
 	double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
 
-	double tmp;
+	double tmp,tmp2;
 
 	int itmp;
 
@@ -206,11 +206,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 				if(xa[rxa]*xa[rxa]+xb[rxb]*xb[rxb]>ya[rya]*ya[rya]+yb[ryb]*yb[ryb]){
 
-	/*				printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
+					printf("+(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
 					printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);
 					printf("{%d%d} %.1f %.1f %.1f",j,k,d[0],d[1],d[2]);
-					printf(" | ");
-*/
+					printf("\n");
+
 
 					vtmp = xa; xa = ya; ya = vtmp;
 					vtmp = xb; xb = yb; yb = vtmp;
@@ -222,12 +222,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					for (i = 0; i<3;++i)
 						d[i] = -d[i];
 
-					/*
 
-					printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
+
+					printf("-(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
 					printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);
 					printf("{%d%d} %.1f %.1f %.1f",j,k,d[0],d[1],d[2]);
-					printf("\n");*/
+					printf("\n");
 				}
 
 
@@ -243,10 +243,23 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],
 					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
 					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
+
+					if(xa[rxa]*xa[rxa]+xb[rxb]*xb[rxb]<=ya[rya]*ya[rya]+yb[ryb]*yb[ryb]){
+
 					tmp
 							= apply0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),
-									fabs(ya[rxa]), fabs(yb[rxb]), d[rxa],
+									fabs(ya[rya]), fabs(yb[ryb]), d[rxa],
 									d[rxb], d[rx]);
+
+					}else{
+
+					tmp
+							= apply0Int(F_par, fabs(ya[rya]), fabs(yb[ryb]),
+									fabs(xa[rxa]), fabs(xb[rxb]), -d[rxa],
+									-d[rxb], -d[rx]);
+					}
+			//		if(fabs(tmp-tmp2)>10e-16)
+			//			printf("wtf");
 					//	 printf("%d%d|%.2f\n",j,k,tmp);
 				} else {
 					//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],
diff --git a/src/save1.mat b/src/save1.mat
new file mode 100644
index 0000000..89b1459
Binary files /dev/null and b/src/save1.mat differ
diff --git a/src/save10.mat b/src/save10.mat
new file mode 100644
index 0000000..d246641
Binary files /dev/null and b/src/save10.mat differ
diff --git a/src/save11.mat b/src/save11.mat
new file mode 100644
index 0000000..093cc3d
Binary files /dev/null and b/src/save11.mat differ
diff --git a/src/save12.mat b/src/save12.mat
new file mode 100644
index 0000000..095aa94
Binary files /dev/null and b/src/save12.mat differ
diff --git a/src/save13.mat b/src/save13.mat
new file mode 100644
index 0000000..e028544
Binary files /dev/null and b/src/save13.mat differ
diff --git a/src/save2.mat b/src/save2.mat
new file mode 100644
index 0000000..42820ed
Binary files /dev/null and b/src/save2.mat differ
diff --git a/src/save3.mat b/src/save3.mat
new file mode 100644
index 0000000..cade101
Binary files /dev/null and b/src/save3.mat differ
diff --git a/src/save4.mat b/src/save4.mat
new file mode 100644
index 0000000..2748262
Binary files /dev/null and b/src/save4.mat differ
diff --git a/src/save5.mat b/src/save5.mat
new file mode 100644
index 0000000..e4d3227
Binary files /dev/null and b/src/save5.mat differ
diff --git a/src/save6.mat b/src/save6.mat
new file mode 100644
index 0000000..1a71948
Binary files /dev/null and b/src/save6.mat differ
diff --git a/src/save7.mat b/src/save7.mat
new file mode 100644
index 0000000..8203440
Binary files /dev/null and b/src/save7.mat differ
diff --git a/src/save8.mat b/src/save8.mat
new file mode 100644
index 0000000..297ebba
Binary files /dev/null and b/src/save8.mat differ
diff --git a/src/save9.mat b/src/save9.mat
new file mode 100644
index 0000000..4ce9f2c
Binary files /dev/null and b/src/save9.mat differ
diff --git a/src/test_solveError_Grid.m b/src/test_solveError_Grid.m
new file mode 100755
index 0000000..4c745b6
--- /dev/null
+++ b/src/test_solveError_Grid.m
@@ -0,0 +1,57 @@
+
+dataIso =[];
+dataAniso =[];
+maxSize = 500;
+Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
+% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren
+% Netz = 'exmpl_2DQuad';
+
+%% Anisotrope Verfeinerung
+disp 'Anisotrop'
+load(Netz)
+
+ref = 0;
+while size(elements,1)<maxSize
+  tic
+  ref = ref+1; 
+
+  %Netz Isotrop verfeinern (4 Teile)
+  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
+  
+  %Matrix mit feinem Netz berechnen
+  A_fine = mex_build_A(coordinates_fine,elements_fine);
+  
+  %rechte Seite Aufbauen
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+  
+  %Lösungsvektor bestimmen
+  x_fine = A_fine\b;
+  
+  %Energienorm^2 bestimmen
+  xe_fine = x_fine'*A_fine*x_fine;
+
+  % Welche Elemente sollen verfeinert werden
+  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
+  
+  % Wie soll verfeinert werden
+  marked = mark(x_fine(f2s)',ind,0.5,0.5);
+
+  %Daten zur Auswertung zwischenspeichern (testweise?)
+  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
+  
+  %Sicherheit falls keine Elemente markiert sind
+  if(sum(marked~=1)==0)
+    disp 'Abgebrochen weil nicht mehr verfeinert wird!'
+    break;
+  end
+  
+  %Verfeinern
+  [coordinates, elements] = refineQuad(coordinates,elements,marked);
+  file = ['save' int2str(ref)];
+  save ([file], 'coordinates', 'elements')
+  toc
+end
+dataAniso
+
+clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine
+
diff --git a/src/test_solveError_Grid2.m b/src/test_solveError_Grid2.m
new file mode 100755
index 0000000..f90d5af
--- /dev/null
+++ b/src/test_solveError_Grid2.m
@@ -0,0 +1,126 @@
+
+dataIso =[];
+dataAniso =[];
+maxSize = 500;
+fileName = './test_save';
+Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
+% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren
+% Netz = 'exmpl_2DQuad';
+
+
+% %% Anisotrope Verfeinerung
+% disp 'Anisotrop1'
+% load(Netz)
+% 
+% ref = 0;
+% while size(elements,1)<maxSize
+%   tic
+%   ref = ref+1; 
+%   
+% 
+% 
+%   %Netz Isotrop verfeinern (4 Teile)
+%   [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
+%   
+%   %Matrix mit feinem Netz berechnen
+%   A_fine = mex_build_A(coordinates_fine,elements_fine);
+%   
+%   %rechte Seite Aufbauen
+%   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+%   
+%   %Lösungsvektor bestimmen
+%   x_fine = A_fine\b;
+%   
+%   %Energienorm^2 bestimmen
+%   xe_fine = x_fine'*A_fine*x_fine;
+% 
+%   % Welche Elemente sollen verfeinert werden
+%   ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
+%   
+%   % Wie soll verfeinert werden
+%   marked = mark(x_fine(f2s)',ind,0.5,0.5);
+% 
+%   %Daten zur Auswertung zwischenspeichern (testweise?)
+%   dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
+%   
+%   %Sicherheit falls keine Elemente markiert sind
+%   if(sum(marked~=1)==0)
+%     disp 'Abgebrochen weil nicht mehr verfeinert wird!'
+%     break;
+%   end
+%   
+%   %Verfeinern
+%   %[coordinates, elements] = refineQuad(coordinates,elements,marked);
+%   file = ['save' int2str(ref)];
+%   load(file);
+%   toc
+% end
+% d1 = dataAniso
+
+
+%% Anisotrope Verfeinerung
+disp 'Anisotrop2'
+load(Netz)
+
+ref = 0;
+while size(elements,1)<maxSize
+  tic
+  ref = ref+1; 
+  
+
+
+  %Netz Isotrop verfeinern (4 Teile)
+  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
+  
+  %Matrix mit feinem Netz berechnen
+  A_fine = mex_build_As1(coordinates_fine,elements_fine);
+  
+  %rechte Seite Aufbauen
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+  
+  %Lösungsvektor bestimmen
+  x_fine = A_fine\b;
+  
+  %Energienorm^2 bestimmen
+  xe_fine = x_fine'*A_fine*x_fine;
+
+  % Welche Elemente sollen verfeinert werden
+  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
+  
+  % Wie soll verfeinert werden
+  marked = mark(x_fine(f2s)',ind,0.5,0.5);
+
+  %Daten zur Auswertung zwischenspeichern (testweise?)
+  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
+  
+  %Sicherheit falls keine Elemente markiert sind
+  if(sum(marked~=1)==0)
+    disp 'Abgebrochen weil nicht mehr verfeinert wird!'
+    break;
+  end
+  
+  %Verfeinern
+  %[coordinates, elements] = refineQuad(coordinates,elements,marked);
+  file = ['save' int2str(ref)];
+  load(file);
+  toc
+end
+d2 = dataAniso
+
+clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine
+
+
+%% Alles Zeichnen
+figure (1)
+loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));
+
+legend('Isotrop','Anisotrop')
+xlabel('Elemente')
+ylabel('Fehler')
+
+figure (2)
+loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));
+
+legend('Isotrop','Anisotrop')
+xlabel('Elemente')
+ylabel('EnergieNorm')