[src] weitere Anpassungen an die test_calc Files

[src] erste Versuche den QuadAchse NAN Fehler zu finden

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

diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp
index 3dbe4c564c30cf9f659f5c6bee837a2dcbd569bc..84d81755fe9e0caab28816b6780ac4333d625b9e 100644 (file)
--- a/src/slpRectangle.cpp
+++ b/src/slpRectangle.cpp
@@ -1,6 +1,7 @@
 #include <iostream>\r
 #include <cmath>\r
 #include <cassert>\r
+#include <limits>\r
 #include <stdlib.h>\r
 \r
 //#include <mex.h>\r
@@ -97,15 +98,20 @@ double G_QY1Y2(double p, double y1, double y2, double x1, double x2, double l) {
 \r
 double G_AY2X2(double y1, double y2, double x1, double x2, double d3) {\r
 \r
-       //       printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
+//     cout << "("     << y1 << "," << y2 << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl;\r
 \r
        double hlp = ((x1 - y1) * (x1 - y1) + (x2 - y2) * (x2 - y2) + d3 * d3);\r
 \r
        double sol = sqrt(hlp);\r
 \r
+\r
        if ((x2 - y2) != 0)\r
                sol += (x2 - y2) * log(sqrt(hlp) - (x2 - y2));\r
 \r
+       if (sol!=sol||sol==numeric_limits<double>::infinity())\r
+               cerr << "G_AY2X2: " << sol << " isn't a Number. ("\r
+                       << y1 << "," << y2 << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl;\r
+\r
        return sol;\r
 }\r
 \r
@@ -274,19 +280,46 @@ double calcParIntQY1X1(double b, double d, double t, double v, double d1,
 \r
        for (int i = 0; i < 5; ++i) {\r
                for (int j = 0; j < 5; ++j) {\r
+\r
+\r
                        sol\r
-                                       += G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i],\r
-                                                       d, d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
+                                       += G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], d, d3)\r
+                                       * t * b * gauss_w5[i] * gauss_w5[j];\r
                        sol\r
-                                       -= G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i],\r
-                                                       0, d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
+                                       -= G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], 0, d3)\r
+                                       * t * b * gauss_w5[i] * gauss_w5[j];\r
                        sol -= G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], d, d3)\r
                                        * t * b * gauss_w5[i] * gauss_w5[j];\r
                        sol += G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], 0, d3)\r
                                        * t * b * gauss_w5[i] * gauss_w5[j];\r
+\r
+                       if (sol!=sol||sol==numeric_limits<double>::infinity()){\r
+                               cout << "->(" << i << "," << j << ")" << endl;\r
+                               cout << "-|(" <<\r
+                                               G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], d, d3)\r
+                                                                                       * t * b * gauss_w5[i] * gauss_w5[j]\r
+                               << "," <<\r
+                                               G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], 0, d3)\r
+                                                                                       * t * b * gauss_w5[i] * gauss_w5[j]\r
+                               << "," <<\r
+                                               G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], d, d3)\r
+                                                                                       * t * b * gauss_w5[i] * gauss_w5[j]\r
+                               << "," <<\r
+                                               G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], 0, d3)\r
+                                                                                       * t * b * gauss_w5[i] * gauss_w5[j]\r
+                               << ")" << endl;\r
+                               cout << "||"<< t * b * gauss_w5[i] * gauss_w5[j] <<"||-----------------" << endl;\r
+                               cout << endl;\r
+                               return sol;\r
+                       }\r
+\r
                }\r
        }\r
 \r
+       if (sol!=sol||sol==numeric_limits<double>::infinity())\r
+               cout << "calcParIntQY1X1: " << sol << " isn't a Number. ("\r
+                       << b << "," << d << ")(" << t << "," << v << ")(" << d1 << "," << d2 << "," << d3 <<")" << endl;\r
+\r
        return sol;\r
 }\r
 \r

diff --git a/src/test_calcInt1.m b/src/test_calcInt1.m
index ad6da4796c99965fc5b3a21471245e308b6121c8..b342959df87adb077fb5f9a412221f7537481a3d 100644 (file)
--- a/src/test_calcInt1.m
+++ b/src/test_calcInt1.m
@@ -1 +1 @@
-\r\rcoo = @(h,diff)[0 0 0;1 0 0; 1 1 0;0 1  0; 0 0 0 ; 0+h 0 0; 0+h h 0 ; 0 h 0]+[zeros(4,3);repmat(diff,4,1)];\relements=[1 2 3 4;5 6 7 8];\rneigh = zeros(2,8);\rdat = [];\r\rdiff = [2 0 0];\r\r%% Laage Übersicht\rfigure(1)\rh = 1;\rcoordinates=coo(h,diff)\rcurrent = coordinates(elements(1,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'g');\rhold on\rcurrent = coordinates(elements(2,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'b');\r\rh = h/2;\rcoordinates=coo(h,diff);\rcurrent = coordinates(elements(2,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'y');\rhold off\rlegend('Element1(h)','Element2(h)','Element2(h/2)');\rtitle('Laage der Elemente');\r\r%% Integrale bei kleiner werdenden Elementen\r\rh = 2;\rfor I = 1:50\r\r  h = h/2;\r  coordinates=coo(h,diff);\r  \r  A0 = mex_build_AU(coordinates,elements,0,0);\r  A2 = mex_build_AU(coordinates,elements,1,2);\r  A1 = mex_build_AU(coordinates,elements,1,1);\r  I\r  dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)];\r    dat(dat<0)=0;\rend\r\rfigure(2)\rloglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4)));\r\rlegend('Analytisch','Quad Element','Element vertauschen','location','southeast');\rxlabel('Elementgroesse (kürzeste Seite)');\rylabel('Integral');\rtitle('Integral bei kleiner werdenden Element');\r\r%% Netzverfeinerung ab einem Wert\r\rh =  (1/2)^25;\rcoordinates=coo(h,diff);\r\rA_loadMesh(coordinates,elements,neigh);\rdatA=[];\rfor I = 1:5\r  datA(I,:) = A_stepIso(5,[0 1 3 2]);\r%   datA(I,:) = A_stepAniso(6,[0 1 3 2],1,0);\r  I\rend\rfigure(3)\rloglog(datA(2:end,1),datA(2:end,3),datA(2:end,1),datA(2:end,6),datA(2:end,1),datA(2:end,12),datA(2:end,1),datA(2:end,9))\r\rlegend('Analytisch','Element vertauschen','Quad Element','Quad Achse','location','northeast');\rxlabel('Anzahl der Elemente');\rylabel('mu Schaetzer');\rtitle('mu Schaetzer mit "schlechtem" Startnetz');\r\rfigure(4)\rloglog(datA(:,1),datA(:,4),datA(:,1),datA(:,7),datA(:,1),datA(:,13),datA(:,1),datA(:,10))\r\rlegend('Analytisch','Element vertauschen','Quad Element','Quad Achse','location','southeast');\rxlabel('Anzahl der Elemente');\rylabel('EnergieNorm ^2 ');\rtitle('EnergieNorm ^2 mit "schlechtem" Startnetz');\r\rdatA\r
\ No newline at end of file
+\r\rcoo = @(h,diff)[0 0 0;1 0 0; 1 1 0;0 1  0; 0 0 0 ; 0+h 0 0; 0+h h 0 ; 0 h 0]+[zeros(4,3);repmat(diff,4,1)];\relements=[1 2 3 4;5 6 7 8];\rneigh = zeros(2,8);\rdat = [];\r\rdiff = [2 0 0];\r\r%% Laage Übersicht\rfigure(1)\rh = 1;\rcoordinates=coo(h,diff)\rcurrent = coordinates(elements(1,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'g');\rhold on\rcurrent = coordinates(elements(2,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'b');\r\rh = h/2;\rcoordinates=coo(h,diff);\rcurrent = coordinates(elements(2,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'y');\rhold off\rlegend('Element1(h)','Element2(h)','Element2(h/2)');\rtitle('Laage der Elemente');\r\r%% Integrale bei kleiner werdenden Elementen\r\rh = 2;\rfor I = 1:50\r\r  h = h/2;\r  coordinates=coo(h,diff);\r  \r  A0 = mex_build_AU(coordinates,elements,0,0);\r  A2 = mex_build_AU(coordinates,elements,1,2);\r  A1 = mex_build_AU(coordinates,elements,1,1);\r  I\r  dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)];\r    dat(dat<0)=0;\rend\r\rfigure(2)\rloglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4)));\r\rlegend('Analytisch','Quad Element','Element vertauschen','location','best');\rxlabel('Elementgroesse (kürzeste Seite)');\rylabel('Integral');\rtitle('Integral bei kleiner werdenden Element');\r\r%% Netzverfeinerung mit schlechtem Netz\r\rh =  (1/2)^25;\rcoordinates=coo(h,diff);\r\rA_loadMesh(coordinates,elements,neigh);\rdatA=[];\rfor I = 1:5\r  datA(I,:) = A_stepIso(1,[0 3 2]);\r%   datA(I,:) = A_stepAniso(1,[0 3 2],1,0);\r  I\rend\rfigure(3)\rloglog(datA(:,1),datA(:,3),datA(:,1),datA(:,9),datA(:,1),datA(:,6))\r\rlegend('Analytisch','Quad Element','Quad Achse','location','best');\rxlabel('Anzahl der Elemente');\rylabel('mu Schaetzer');\rtitle('mu Schaetzer mit "schlechtem" Startnetz');\r\rfigure(4)\rloglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7))\r\rlegend('Analytisch','Quad Element','Quad Achse','location','best');\rxlabel('Anzahl der Elemente');\rylabel('EnergieNorm ^2 ');\rtitle('EnergieNorm ^2 mit "schlechtem" Startnetz');\r\rdatA\r
\ No newline at end of file

diff --git a/src/test_calcInt2.m b/src/test_calcInt2.m
index c86ccfcf1878963c4f4e5af7fe977d892fe0e3e7..34a9e2d9d2d06bfd24768727c74510622636a128 100644 (file)
--- a/src/test_calcInt2.m
+++ b/src/test_calcInt2.m
@@ -1,15 +1,16 @@
 
 
+coo = @(h,diff)[0 0 0;1 0 0; 1 1 0;0 1  0; 0 0 0 ; 0+h 0 0; 0+h 1 0 ; 0 1 0]+[zeros(4,3);repmat(diff,4,1)];
 elements=[1 2 3 4;5 6 7 8];
 neigh = zeros(2,8);
-
 dat = [];
 
-diff = 10;
+diff = [5 0 0];
 
-h = 1/2;
-coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; diff 0 0; diff+h 0 0; diff+h 1 0 ; diff 1 0];
+%% Laage Übersicht
 figure(1)
+h = 1;
+coordinates=coo(h,diff)
 current = coordinates(elements(1,[1:4,1])',:);
 fill3(current(:,1),current(:,2),current(:,3),'g');
 hold on
@@ -17,30 +18,63 @@ current = coordinates(elements(2,[1:4,1])',:);
 fill3(current(:,1),current(:,2),current(:,3),'b');
 
 h = h/2;
-coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; diff 0 0; diff+h 0 0; diff+h 1 0 ; diff 1 0];
+coordinates=coo(h,diff);
 current = coordinates(elements(2,[1:4,1])',:);
 fill3(current(:,1),current(:,2),current(:,3),'y');
 hold off
-legend('Element1(h/2)','Element2(h/2)','Element2(h/4)')
+legend('Element1(h)','Element2(h)','Element2(h/2)');
+title('Laage der Elemente');
+
+%% Integrale bei kleiner werdenden Elementen
 
 h = 2;
 for I = 1:50
 
   h = h/2;
-coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; diff 0 0; diff+h 0 0; diff+h 1 0 ; diff 1 0];
+  coordinates=coo(h,diff);
   
   A0 = mex_build_AU(coordinates,elements,0,0);
   A2 = mex_build_AU(coordinates,elements,1,2);
-  A3 = mex_build_AU(coordinates,elements,1,3);
+  A1 = mex_build_AU(coordinates,elements,1,3);
   I
-  dat(I,1:4) = [h A0(1,2) A2(1,2) A3(1,2)];
-%    dat(dat<0)=0;
+  dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)];
+    dat(dat<0)=0;
 end
 
 figure(2)
 loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4)));
 
-legend('Analytisch','Quad Element','Quad Achse');
+legend('Analytisch','Quad Element','Quad Achse','location','best');
+xlabel('Elementgroesse (kürzeste Seite)');
+ylabel('Integral');
+title('Integral bei kleiner werdenden Element');
+
+%% Netzverfeinerung mit schlechtem Netz
+
+h =  (1/2)^25;
+coordinates=coo(h,diff);
+
+A_loadMesh(coordinates,elements,neigh);
+datA=[];
+for I = 1:5
+  datA(I,:) = A_stepIso(1,[0 3 2]);
+%   datA(I,:) = A_stepAniso(1,[0 1 3 2],1,0);
+  I
+end
+figure(3)
+loglog(datA(:,1),datA(:,3),datA(:,1),datA(:,9),datA(:,1),datA(:,6))
+
+legend('Analytisch','Quad Element','Quad Achse','location','best');
+xlabel('Anzahl der Elemente');
+ylabel('mu Schaetzer');
+title('mu Schaetzer mit "schlechtem" Startnetz');
+
+figure(4)
+loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7))
+
+legend('Analytisch','Quad Element','Quad Achse','location','best');
+xlabel('Anzahl der Elemente');
+ylabel('EnergieNorm ^2 ');
+title('EnergieNorm ^2 mit "schlechtem" Startnetz');
 
-xlabel('Elementgroesse');
-ylabel('Integral');
\ No newline at end of file
+datA

diff --git a/src/test_calcInt3.m b/src/test_calcInt3.m
index f83514e215f913436c4c64ea791dc5d7aac4ae12..d530636e6c34eb1ef5dbb5d2965a6d574a034849 100644 (file)
--- a/src/test_calcInt3.m
+++ b/src/test_calcInt3.m
@@ -1,13 +1,16 @@
 
 
+coo = @(h,diff)[0 0 0;h 0 0; h 1 0;0 1  0; 0 0 0 ; 0+h 0 0; 0+h 1 0 ; 0 1 0]+[zeros(4,3);repmat(diff,4,1)];
 elements=[1 2 3 4;5 6 7 8];
 neigh = zeros(2,8);
-
 dat = [];
 
-h = 1/2;
-coordinates=[0 0 0;h 0 0; h 1 0;0 1  0; 2 0 0; 2+h 0 0; 2+h 1 0 ; 2 1 0];
+diff = [5 0 0];
+
+%% Laage Übersicht
 figure(1)
+h = 1;
+coordinates=coo(h,diff)
 current = coordinates(elements(1,[1:4,1])',:);
 fill3(current(:,1),current(:,2),current(:,3),'g');
 hold on
@@ -15,33 +18,65 @@ current = coordinates(elements(2,[1:4,1])',:);
 fill3(current(:,1),current(:,2),current(:,3),'b');
 
 h = h/2;
-coordinates=[0 0 0;h 0 0; h 1 0;0 1  0; 2 0 0; 2+h 0 0; 2+h 1 0 ; 2 1 0];
+coordinates=coo(h,diff);
 current = coordinates(elements(1,[1:4,1])',:);
-fill3(current(:,1),current(:,2),current(:,3),'y');
-hold on
-current = coordinates(elements(2,[1:4,1])',:);
 fill3(current(:,1),current(:,2),current(:,3),'r');
+current = coordinates(elements(2,[1:4,1])',:);
+fill3(current(:,1),current(:,2),current(:,3),'y');
 hold off
-legend('Element1(h/2)','Element2(h/2)','Element1(h/4)','Element2(h/4)')
+legend('Element1(h)','Element2(h)','Element1(h/2)','Element2(h/2)');
+title('Laage der Elemente');
+
+%% Integrale bei kleiner werdenden Elementen
 
 h = 2;
 for I = 1:50
 
   h = h/2;
-coordinates=[0 0 0;h 0 0; h 1 0;0 1  0; 2 0 0; 2+h 0 0; 2+h 1 0 ; 2 1 0];
+  coordinates=coo(h,diff);
   
   A0 = mex_build_AU(coordinates,elements,0,0);
   A2 = mex_build_AU(coordinates,elements,1,2);
-  A3 = mex_build_AU(coordinates,elements,1,3);
+  A1 = mex_build_AU(coordinates,elements,1,3);
   I
-  dat(I,1:4) = [h A0(1,2) A2(1,2) A3(1,2)];
+  dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)];
     dat(dat<0)=0;
 end
 
 figure(2)
 loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4)));
 
-legend('Analytisch','Quad Element','Quad Achse');
+legend('Analytisch','Quad Element','Quad Achse','location','best');
+xlabel('Elementgroesse (kürzeste Seite)');
+ylabel('Integral');
+title('Integral bei kleiner werdenden Element');
+
+%% Netzverfeinerung mit schlechtem Netz
+
+h =  (1/2)^25;
+coordinates=coo(h,diff);
+
+A_loadMesh(coordinates,elements,neigh);
+datA=[];
+for I = 1:5
+  datA(I,:) = A_stepIso(1,[0 3 2]);
+%   datA(I,:) = A_stepAniso(1,[0 3 2],1,0);
+  I
+end
+figure(3)
+loglog(datA(:,1),datA(:,3),datA(:,1),datA(:,9),datA(:,1),datA(:,6))
+
+legend('Analytisch','Quad Element','Quad Achse','location','best');
+xlabel('Anzahl der Elemente');
+ylabel('mu Schaetzer');
+title('mu Schaetzer mit "schlechtem" Startnetz');
+
+figure(4)
+loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7))
+
+legend('Analytisch','Quad Element','Quad Achse','location','best');
+xlabel('Anzahl der Elemente');
+ylabel('EnergieNorm ^2 ');
+title('EnergieNorm ^2 mit "schlechtem" Startnetz');
 
-xlabel('Elementgroesse');
-ylabel('Integral');
\ No newline at end of file
+datA