[src] neue Test's mit einem Element

[src] mexTest -> läuft MEX überhaupt

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

diff --git a/src/mex_build_AU.cpp b/src/mex_build_AU.cpp
index a57c1972697d18c2646753a106c980ca9574d441..3dc3aae00fd1c2996024b1d8604865d11e0f6fc0 100644 (file)
--- a/src/mex_build_AU.cpp
+++ b/src/mex_build_AU.cpp
@@ -21,14 +21,14 @@
 using namespace std;\r
 \r
 int dimOfVec(double* vec) {\r
-       if (fabs(vec[2]) > M_EPS)\r
+       if (vec[2] != 0)\r
                return 2;\r
-       else if (fabs(vec[1]) > M_EPS)\r
+       else if (vec[1] !=0)\r
                return 1;\r
-       else if (fabs(vec[0]) > M_EPS)\r
+       else if (vec[0] !=0)\r
                return 0;\r
        else {\r
-               cerr << "dimOfVec : (" << vec[0] << vec[1] << vec[2]\r
+               cerr << "dimOfVec : (" << vec[0] << " " << vec[1]<< " " << vec[2]\r
                                << ") alle Werte 0 \n";\r
                return -1;\r
        }\r
@@ -105,22 +105,26 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
        case 0:\r
                ctypeP = calcParIntO0;\r
                ctypeO = calcOrtIntO0;\r
+               break;\r
        case 1:\r
                ctypeP = calcParIntO1;\r
                ctypeO = calcOrtIntO1;\r
+               break;\r
        case 2:\r
                ctypeP = calcParIntO2;\r
                ctypeO = calcOrtIntO2;\r
+               break;\r
        case 3:\r
                ctypeP = calcParIntO3;\r
                ctypeO = calcOrtIntO3;\r
+               break;\r
 \r
        }\r
 \r
        //LageInformationen\r
        int rx, rxa, rxb, ry, rya, ryb;\r
 \r
-       cout << "  Progress: #";\r
+//     cout << "  Progress: #";\r
        count = 0;\r
        //Ausrechnen\r
        for (j = 0; j < em; ++j) {\r
@@ -263,11 +267,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                        A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
                        if(k!=j)\r
                                A[(j * em) + k] = 1. / (4 * M_PI) * tmp;\r
-                       if(count++ > ((em*(em+1))/2)/10){\r
+/*                     if(count++ > ((em*(em+1))/2)/10){\r
                                count = 0;\r
                                cout << "#";\r
                                cout.flush();\r
-                       }\r
+                       }*/\r
 \r
                }\r
 \r

diff --git a/src/mex_hello.cpp b/src/mex_hello.cpp
new file mode 100644 (file)
index 0000000..7911037
--- /dev/null
+++ b/src/mex_hello.cpp
@@ -0,0 +1,9 @@
+#include <iostream>
+#include "mex.h"
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+       std::cout << "Hello" << std::endl;
+
+
+
+}

diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp
index 6a294955c7172c15325fe0229cea0f4b5db461b0..d984cf82e588a4d4e27fc4d26b7a9f5f14cc9790 100644 (file)
--- a/src/slpRectangle.cpp
+++ b/src/slpRectangle.cpp
@@ -283,7 +283,6 @@ double calcParIntQY1X1(double b, double d, double t, double v, double d1,
        }\r
 \r
        return sol;\r
-\r
 }\r
 \r
 double calcParIntQY1(double b, double d, double t, double v, double d1,\r
@@ -299,9 +298,25 @@ double calcOrtIntA(double b, double d, double t, double v, double d1,
 \r
 }\r
 \r
-double calcOrtIntQY1Y2(double b, double d, double t, double v, double d1,\r
+double calcOrtIntQX1X2(double b, double d, double t, double v, double d1,\r
                double d2, double d3) {\r
-       return 0;\r
+       //Fall 2\r
+       double sol = 0;\r
+\r
+       for (int i = 0; i < 5; ++i) {\r
+               for (int j = 0; j < 5; ++j) {\r
+                       sol += G_AY1Y2(-0.5, t + d2, v + d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1)\r
+                                       * b * d * gauss_w5[i] * gauss_w5[j];\r
+                       sol -= G_AY1Y2(-0.5, d2, v + d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1)\r
+                                       * b * d * gauss_w5[i] * gauss_w5[j];\r
+                       sol -= G_AY1Y2(-0.5, t + d2, d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1)\r
+                                       * b * d * gauss_w5[i] * gauss_w5[j];\r
+                       sol += G_AY1Y2(-0.5, d2, d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1)\r
+                                       * b * d * gauss_w5[i] * gauss_w5[j];\r
+               }\r
+       }\r
+\r
+       return sol;\r
 \r
 }\r
 \r
@@ -344,11 +359,11 @@ double calcOrtIntO1(double b, double d, double t, double v, double d1,
                double tmp = 0;\r
 \r
                tmp = b;\r
-               b = t;\r
-               t = tmp;\r
-               tmp = d;\r
-               d = v;\r
+               b = v;\r
                v = tmp;\r
+               tmp = d;\r
+               d = t;\r
+               t = tmp;\r
 \r
                d1 = -d1;\r
                d2 = -d2;\r
@@ -389,15 +404,15 @@ double calcOrtIntO2(double b, double d, double t, double v, double d1,
                double d2, double d3, double eta) {\r
 \r
        //Elmente anordnen\r
-       if (b * b + d * d > t * t + v * v) {\r
+       if (d > t) {\r
                double tmp = 0;\r
 \r
                tmp = b;\r
-               b = t;\r
-               t = tmp;\r
-               tmp = d;\r
-               d = v;\r
+               b = v;\r
                v = tmp;\r
+               tmp = d;\r
+               d = t;\r
+               t = tmp;\r
 \r
                d1 = -d1;\r
                d2 = -d2;\r
@@ -405,7 +420,7 @@ double calcOrtIntO2(double b, double d, double t, double v, double d1,
        }\r
 \r
        if ((b * b + d * d) * eta <= d1 * d1 + d2 * d2 + d3 * d3) {\r
-               return 0;\r
+               return calcOrtIntQX1X2(b, d, t, v, d1, d2, d3);\r
        } else {\r
                return calcOrtIntA(b, d, t, v, d1, d2, d3);\r
        }\r
@@ -442,15 +457,15 @@ double calcOrtIntO3(double b, double d, double t, double v, double d1,
                double d2, double d3, double eta) {\r
 \r
        //Elmente anordnen\r
-       if (b * b + d * d > t * t + v * v) {\r
+       if (d > t) {\r
                double tmp = 0;\r
 \r
                tmp = b;\r
-               b = t;\r
-               t = tmp;\r
-               tmp = d;\r
-               d = v;\r
+               b = v;\r
                v = tmp;\r
+               tmp = d;\r
+               d = t;\r
+               t = tmp;\r
 \r
                d1 = -d1;\r
                d2 = -d2;\r

diff --git a/src/test_calcInt1.m b/src/test_calcInt1.m
new file mode 100644 (file)
index 0000000..5a9e679
--- /dev/null
+++ b/src/test_calcInt1.m
@@ -0,0 +1 @@
+\r\relements=[1 2 3 4;5 6 7 8];\rneigh = zeros(2,8);\r\rdat = [];\r\rh = 1;\rcoordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 2 0 0; 2+h 0 0; 2+h h 0 ; 2 h 0];\rfigure(1)\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=[0 0 0;1 0 0; 1 1 0;0 1  0; 2 0 0; 2+h 0 0; 2+h h 0 ; 2 h 0];\rcurrent = coordinates(elements(2,[1:4,1])',:);\rfill3(current(:,1),current(:,2),current(:,3),'y');\rhold off\rlegend('Element1(h/2)','Element2(h/2)','Element2(h/4)')\r\rh = 2;\rfor I = 1:50\r\r  h = h/2;\rcoordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 2 0 0; 2+h 0 0; 2+h h 0 ; 2 h 0];\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');\r\rxlabel('Elementgroesse');\rylabel('Integral');
\ No newline at end of file

diff --git a/src/test_calcInt2.m b/src/test_calcInt2.m
new file mode 100644 (file)
index 0000000..9320aec
--- /dev/null
+++ b/src/test_calcInt2.m
@@ -0,0 +1,44 @@
+
+
+elements=[1 2 3 4;5 6 7 8];
+neigh = zeros(2,8);
+
+dat = [];
+
+h = 1/2;
+coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 4 0 0; 4+h 0 0; 4+h 1 0 ; 4 1 0];
+figure(1)
+current = coordinates(elements(1,[1:4,1])',:);
+fill3(current(:,1),current(:,2),current(:,3),'g');
+hold on
+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; 4 0 0; 4+h 0 0; 4+h 1 0 ; 4 1 0];
+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)')
+
+h = 2;
+for I = 1:50
+
+  h = h/2;
+coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 4 0 0; 4+h 0 0; 4+h 1 0 ; 4 1 0];
+  
+  A0 = mex_build_AU(coordinates,elements,0,0);
+  A2 = mex_build_AU(coordinates,elements,1,2);
+  A3 = 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;
+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');
+
+xlabel('Elementgroesse');
+ylabel('Integral');
\ No newline at end of file

diff --git a/src/test_calcInt3.m b/src/test_calcInt3.m
new file mode 100644 (file)
index 0000000..f83514e
--- /dev/null
+++ b/src/test_calcInt3.m
@@ -0,0 +1,47 @@
+
+
+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];
+figure(1)
+current = coordinates(elements(1,[1:4,1])',:);
+fill3(current(:,1),current(:,2),current(:,3),'g');
+hold on
+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];
+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');
+hold off
+legend('Element1(h/2)','Element2(h/2)','Element1(h/4)','Element2(h/4)')
+
+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];
+  
+  A0 = mex_build_AU(coordinates,elements,0,0);
+  A2 = mex_build_AU(coordinates,elements,1,2);
+  A3 = 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;
+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');
+
+xlabel('Elementgroesse');
+ylabel('Integral');
\ No newline at end of file

diff --git a/src/test_calcInt4.m b/src/test_calcInt4.m
new file mode 100644 (file)
index 0000000..c2b60d7
--- /dev/null
+++ b/src/test_calcInt4.m
@@ -0,0 +1,46 @@
+
+
+elements=[1 2 3 4;5 6 7 8];
+neigh = zeros(2,8);
+
+dat = [];
+
+h = 1/2;
+coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 2 0 0; 2 0 h; 2 h h ; 2 h 0];
+figure(1)
+current = coordinates(elements(1,[1:4,1])',:);
+fill3(current(:,1),current(:,2),current(:,3),'g');
+hold on
+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; 2 0 0; 2 0 h; 2 h h ; 2 h 0];
+current = coordinates(elements(2,[1:4,1])',:);
+fill3(current(:,1),current(:,2),current(:,3),'r');
+hold off
+legend('Element1(h/2)','Element2(h/2)','Element1(h/4)','Element2(h/4)')
+
+h = 2;
+for I = 1:50
+
+  h = h/2;
+coordinates=[0 0 0;1 0 0; 1 1 0;0 1  0; 2 0 0; 2 0 h; 2 h h ; 2 h 0];
+  
+  A0 = mex_build_AU(coordinates,elements,0,0);
+  A1 = mex_build_AU(coordinates,elements,1,1);
+  A2 = mex_build_AU(coordinates,elements,1,2);
+  I
+  dat(I,1:4) = [h A0(1,2) A1(1,2) A2(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','Element vertauschen','Quad Element','location','southeast');
+
+xlabel('Elementgroesse');
+ylabel('Integral');
\ No newline at end of file