[src] tests für g0 G00 und F_par hinzugefügt

[src] kleine Fehler behoben...

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

diff --git a/src/SLPrecangle.cpp b/src/SLPrecangle.cpp
index 5a414b91390ff3efae3d92e571a75f9dc4c6f12f..d5c9e111acbb1bf55a43cf2481a907861ee76e43 100644 (file)
--- a/src/SLPrecangle.cpp
+++ b/src/SLPrecangle.cpp
@@ -12,13 +12,8 @@
 using namespace std;\r
 \r
 int sign(double);\r
-double arsinh(double);\r
+//double arsinh(double);\r
 \r
-// sol = quadInt((F_par/F_ort),s1,s2,k1,k2,t1,t2,l1,l2,d1,d2,d3);\r
-double quadInt(\r
-               double(*)(double, double, double, double, double, double, double),\r
-               double, double, double, double, double, double, double, double, double,\r
-               double, double);\r
 /*\r
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
        //sicherheitsabfragen zu Datengroessen\r
@@ -166,11 +161,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 int inline sign(double x) {\r
        return x > 0 ? 1 : (x < 0 ? -1 : 0);\r
 }\r
-\r
+/*\r
 double inline arsinh(double x) {\r
        return log(x + sqrt(x * x + 1));\r
 }\r
-\r
+*/\r
 //y-x muss != 0 sein\r
 double g0(double p, double y, double x, double l) {\r
        //printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);\r
@@ -187,7 +182,7 @@ double g0(double p, double y, double x, double l) {
                else if (p == -0.5)\r
                        sol = asinh((y - x) / fabs(l));\r
                else if (p == -1)\r
-                       sol = atan((y - x) / fabs(l));\r
+                       sol = atan((y - x) / fabs(l))/2;\r
                else if (p == -1.5)\r
                        sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l));\r
                else\r
@@ -229,7 +224,9 @@ double G00(double p, double y1, double y2, double x1, double x2, double l) {
                                        sqrt((y2 - x2) * (y2 - x2) + l * l));\r
                sol /= 2 * p + 2;\r
        } else {\r
-               mexErrMsgTxt("no case for p defined");\r
+               sol = NAN;\r
+               printf("warning in G00: no case for p defined\n");\r
+               //mexErrMsgTxt("no case for p defined\n");\r
        }\r
 \r
        return sol;\r
@@ -242,14 +239,14 @@ double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
        double sol = (x1 - y1 - d1) * (x2 - y2 - d2);\r
 \r
        if (sol != 0)\r
-               sol *= slpADLO(x1, x2, y1 + d1, y2 + d2, d3);\r
+               sol *= G00(-0.5,x1, x2, y1 + d1, y2 + d2, d3);\r
 \r
        if ((x1 - y1 - d1) != 0)\r
-               sol -= (x1 - y1 - d1) * compute_g0(0.5, x1, y1 + d1,\r
+               sol -= (x1 - y1 - d1) * g0(0.5, x1, y1 + d1,\r
                                sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));\r
 \r
        if ((x2 - y2 - d2) != 0)\r
-               sol -= (x2 - y2 - d2) * compute_g0(0.5, x2, y2 + d2,\r
+               sol -= (x2 - y2 - d2) * g0(0.5, x2, y2 + d2,\r
                                sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));\r
 \r
        double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
@@ -258,11 +255,16 @@ double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
        return sol;\r
 }\r
 \r
-double inline quadInt(\r
+double quadInt(\r
                double(*f)(double, double, double, double, double, double, double),\r
-               double s1, double s2, double k1, double k2, double t1, double t2,\r
-               double l1, double l2, double d1, double d2, double d3) {\r
+               double a, double b, double c, double d, double r, double t,\r
+               double u, double v, double d1, double d2, double d3) {\r
 \r
+   return f(b,d,t,v,d1,d2,d3)-f(b,d,t,u,d1,d2,d3)-f(b,d,r,v,d1,d2,d3)+f(b,d,r,u,d1,d2,d3)\r
+    -f(b,c,t,v,d1,d2,d3)+f(b,c,t,u,d1,d2,d3)+f(b,c,r,v,d1,d2,d3)-f(b,c,r,u,d1,d2,d3)\r
+    -f(a,d,t,v,d1,d2,d3)+f(a,d,t,u,d1,d2,d3)+f(a,d,r,v,d1,d2,d3)-f(a,d,r,u,d1,d2,d3)\r
+    +f(a,c,t,v,d1,d2,d3)-f(a,c,t,u,d1,d2,d3)-f(a,c,r,v,d1,d2,d3)+f(a,c,r,u,d1,d2,d3);\r
+/*\r
        return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f(\r
                        k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1,\r
                        s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2,\r
@@ -271,8 +273,14 @@ double inline quadInt(
                        d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1,\r
                        d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2,\r
                        d3) + f(s1, s2, t1, t2, d1, d2, d3);\r
+\r
+                       */\r
 }\r
 \r
+\r
+\r
+\r
+\r
 double slpADLO(double y1, double y2, double x1, double x2, double a) {\r
        double G3 = 0;\r
        double gL = 0;\r

diff --git a/src/SLPrecangle.hpp b/src/SLPrecangle.hpp
index 6f801507e1cf8962374354c4434e430433a54707..0c915487957b94a8bcc4af575af596656c2408fe 100644 (file)
--- a/src/SLPrecangle.hpp
+++ b/src/SLPrecangle.hpp
@@ -18,4 +18,10 @@ double slpADLO(double, double, double, double, double);
 double compute_g0(double, double, double, double);
 double FLO_plane(double, double, double, double, double, double, double);
 
+// sol = quadInt((F_par/F_ort),s1,s2,k1,k2,t1,t2,l1,l2,d1,d2,d3);
+double quadInt(
+               double(*)(double, double, double, double, double, double, double),
+               double, double, double, double, double, double, double, double, double,
+               double, double);
+
 #endif

diff --git a/src/build_A.cpp b/src/build_A.cpp
index f78998570fb93de2cea1cb053ad342c630694c77..03c5054f3d293f784d8138412ca2b64171f8849e 100644 (file)
--- a/src/build_A.cpp
+++ b/src/build_A.cpp
@@ -4,35 +4,13 @@
 #include "mex.h"\r
 #include <stdlib.h>\r
 \r
+#include "SLPrecangle.hpp"\r
+\r
 #define my_PI 3.141592653589793\r
 #define EPS 0.02\r
 \r
 using namespace std;\r
 \r
-int sign(double);\r
-double arsinh(double);\r
-\r
-// sol = g0(p,y,x,l);\r
-double g0(double, double, double, double);\r
-// sol = G00(p,y1,y2,x1,x2,l);\r
-double G00(double, double, double, double, double, double);\r
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3);\r
-double F_par(double, double, double, double, double, double, double);\r
-// sol = F_ort(x1,x2,y1,y2,d1,d2,d3);\r
-//double F_ort(double, double, double, double, double, double, double);\r
-\r
-\r
-double slpADLO(double, double, double, double, double);\r
-double compute_g0(double, double, double, double);\r
-double FLO_plane(double, double, double, double, double, double, double);\r
-\r
-\r
-// sol = quadInt((F_par/F_ort),s1,s2,k1,k2,t1,t2,l1,l2,d1,d2,d3);\r
-double quadInt(\r
-               double(*)(double, double, double, double, double, double, double),\r
-               double, double, double, double, double, double, double, double, double,\r
-               double, double);\r
-\r
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
        //sicherheitsabfragen zu Datengroessen\r
        if (nrhs != 2)\r
@@ -138,8 +116,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                                if (rx = 2) {\r
                                        //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x1[0],x1[1],x2[0],x3[1],y1[0],y1[1],y2[0],y3[1],d[0],d[1],d[2]);\r
                                        //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
-                                       tmp = quadInt(FLO_plane, x1[0], x1[1], x2[0], x3[1], y1[0],\r
-                                                       y1[1], y2[0], y3[1], d[0], d[1], d[2]);\r
+                                       tmp = quadInt(F_par, x1[0], x2[0], x1[1], x3[1], y1[0],\r
+                                                       y2[0], y1[1], y3[1], d[0], d[1], d[2]);\r
                                        //       printf("%d%d|%.2f\n",j,k,tmp);\r
                                        A[(k * em) + j] = 1. / (4 * my_PI) * tmp;\r
 \r
@@ -175,198 +153,3 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 \r
        return;\r
 }\r
-\r
-int inline sign(double x) {\r
-       return x > 0 ? 1 : (x < 0 ? -1 : 0);\r
-}\r
-\r
-double inline arsinh(double x) {\r
-       return log(x + sqrt(x * x + 1));\r
-}\r
-\r
-//y-x muss != 0 sein\r
-double g0(double p, double y, double x, double l) {\r
-       //printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);\r
-\r
-       double sol = 0;\r
-\r
-       if (l != 0) {\r
-               if (p == 0.5) {\r
-                       sol = (y - x) / 2 * sqrt((y - x) * (y - x) + l * l) + l * l / 2\r
-                                       * arsinh((y - x) / fabs(l));\r
-                       //          printf("%.2f |",sol);\r
-               } else if (p == 0)\r
-                       sol = y - x;\r
-               else if (p == -0.5)\r
-                       sol = asinh((y - x) / fabs(l));\r
-               else if (p == -1)\r
-                       sol = atan((y - x) / fabs(l));\r
-               else if (p == -1.5)\r
-                       sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l));\r
-               else\r
-                       sol = (y - x) * pow((y - x) * (y - x) + l * l, p) + 2 * p * l * l\r
-                                       * g0(p - 1, y, x, l) / (2 * p + 1);\r
-       } else {\r
-               if (p == -0.5)\r
-                       sol = sign(y - x) * log(fabs(y - x));\r
-               else\r
-                       sol = (y - x) * pow(fabs(y - x), 2 * p) / (2 * p + 1);\r
-       }\r
-\r
-       return sol;\r
-}\r
-\r
-double G00(double p, double y1, double y2, double x1, double x2, double l) {\r
-       //      printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);\r
-\r
-       double sol = 0;\r
-       if (p == -1.5) {\r
-               if (l == 0) {\r
-                       sol = -sqrt((y1 - x1) * (y1 - x1) + (y2 - x2) * (y2 - x2)) / ((y1\r
-                                       - x1) * (y2 - x2));\r
-               } else {\r
-                       sol = sign((y1 - x1) * (y2 - x2)) / (2 * fabs(l)) * acos(\r
-                                       -2 * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2) / (((y1\r
-                                                       - x1) * (y1 - x1) + l * l) * ((y2 - x2) * (y2 - x2)\r
-                                                       + l * l)) + 1);\r
-               }\r
-\r
-       } else if (p == -0.5) {\r
-               if (l != 0)\r
-                       sol = 2 * p * l * l * G00(p - 1, y1, y2, x1, x2, l);\r
-               if ((y1 - x1) != 0)\r
-                       sol += (y1 - x1) * g0(p, y2, x2,\r
-                                       sqrt((y1 - x1) * (y1 - x1) + l * l));\r
-               if ((y2 - x2) != 0)\r
-                       sol += (y2 - x2) * g0(p, y1, x1,\r
-                                       sqrt((y2 - x2) * (y2 - x2) + l * l));\r
-               sol /= 2 * p + 2;\r
-       } else {\r
-               mexErrMsgTxt("no case for p defined");\r
-       }\r
-\r
-       return sol;\r
-}\r
-\r
-double F_par(double x1, double x2, double y1, double y2, double d1, double d2,\r
-               double d3) {\r
-\r
-       //       printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
-       double sol = (x1 - y1 - d1) * (x2 - y2 - d2);\r
-\r
-       if (sol != 0)\r
-               sol *= slpADLO( x1, x2, y1 + d1, y2 + d2, d3);\r
-\r
-       if ((x1 - y1 - d1) != 0)\r
-               sol -= (x1 - y1 - d1) * compute_g0(0.5, x1, y1 + d1,\r
-                               sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));\r
-\r
-       if ((x2 - y2 - d2) != 0)\r
-               sol -= (x2 - y2 - d2) * compute_g0(0.5, x2, y2 + d2,\r
-                               sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));\r
-\r
-       double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
-                       - d2) + d3 * d3);\r
-       sol += 1. / 3 * hlp * sqrt(hlp);\r
-       return sol;\r
-}\r
-\r
-double inline quadInt(\r
-               double(*f)(double, double, double, double, double, double, double),\r
-               double s1, double s2, double k1, double k2, double t1, double t2,\r
-               double l1, double l2, double d1, double d2, double d3) {\r
-\r
-       return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f(\r
-                       k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1,\r
-                       s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2,\r
-                       t1, l2, d1, d2, d3) - f(k1, s2, t1, t2, d1, d2, d3) - f(s1, k2, l1,\r
-                       l2, d1, d2, d3) + f(s1, k2, l1, t2, d1, d2, d3) + f(s1, k2, t1, l2,\r
-                       d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1,\r
-                       d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2,\r
-                       d3) + f(s1, s2, t1, t2, d1, d2, d3);\r
-}\r
-\r
-double slpADLO(double y1, double y2, double x1, double x2, double a) {\r
-       double G3 = 0;\r
-       double gL = 0;\r
-       double gK = 0;\r
-       double tmp;\r
-\r
-       tmp = y1 - x1;\r
-       if (fabs(tmp) >= EPS * y1) {\r
-               tmp = sqrt(y1 * y1 + x1 * x1 + a * a - 2 * y1 * x1);\r
-               gL = compute_g0(-0.5, y2, x2, tmp);\r
-       }\r
-       tmp = y2 - x2;\r
-       if (fabs(tmp) >= EPS * y2) {\r
-               tmp = sqrt(y2 * y2 + x2 * x2 + a * a - 2 * y2 * x2);\r
-               gK = compute_g0(-0.5, y1, x1, tmp);\r
-       }\r
-       if (fabs(a * a) > EPS) {\r
-               if ((y1 - x1) * (y2 - x2) * a >= 0)\r
-                       tmp = 1.;\r
-               else\r
-                       tmp = -1.;\r
-\r
-               G3 = tmp * acos(\r
-                               (-2. * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2)) / (((y1\r
-                                               - x1) * (y1 - x1) + a * a) * ((y2 - x2) * (y2 - x2) + a\r
-                                               * a)) + 1.) / (2. * a);\r
-       }\r
-\r
-       return (y1 - x1) * gL + (y2 - x2) * gK - a * a * G3;\r
-}\r
-\r
-double compute_g0(double p, double y, double x, double a) {\r
-       int sp = (int) 2 * (p - EPS); // MK (p-EPS) instead of (p+EPS)\r
-       //printf("\n compute_g0, p = %lf, sp = %d\n",p,sp);\r
-       assert(\r
-                       p == 0 || (p == -0.5) || (p == -1) || (p == -1.5) || (fabs(a)\r
-                                       <= EPS));\r
-       if (fabs(a) <= EPS) {\r
-               // printf("\n a < eps\n");\r
-               switch (sp) {\r
-               case 0:\r
-                       return y - x;\r
-               case -1:\r
-                       return log(fabs(x - y)) * (y - x) / fabs(y - x);\r
-               case -2:\r
-                       return -(y - x) / fabs(y - x) / fabs(y - x);\r
-               case -3:\r
-                       return -0.5 * (y - x) / fabs(y - x) / fabs(y - x) / fabs(y - x);\r
-               }\r
-       } else {\r
-               //  printf("\n a > eps\n");\r
-               switch (sp) {\r
-               case 0:\r
-                       return y - x;\r
-               case -1:\r
-                       return asinh((y - x) / fabs(a));\r
-               case -2:\r
-                       return atan((y - x) / fabs(a));\r
-               case -3:\r
-                       return (y - x) * pow((x * x + y * y + a * a - 2 * x * y), -0.5)\r
-                                       / (a * a);\r
-               default:\r
-                       printf("p must be either 0, -1/2, -1 or -3/2.");\r
-                       return NAN;\r
-               }\r
-       }\r
-}\r
-\r
-double FLO_plane(double x1, double x2, double y1, double y2, double delta1, double delta2, double a){\r
-  double yd1  = y1+delta1;\r
-  double yd2  = y2+delta2;\r
-  double tmp1 = x1-y1-delta1;\r
-  double tmp2 = x2-y2-delta2;\r
-  double tmp3 = sqrt( tmp2*tmp2 + a*a );\r
-  double tmp4 = sqrt( tmp1*tmp1 + a*a );\r
-  double tmp5 = pow(tmp1*tmp1 + tmp2*tmp2 + a*a,3./2.);\r
-  double rval = 0;\r
-\r
-  rval = tmp1*tmp2*slpADLO(x1,x2,yd1,yd2,a)\r
-    - tmp1*compute_g0(-0.5,x1,yd1,tmp3)\r
-    - tmp2*compute_g0(-0.5,x2,yd2,tmp4)\r
-    + tmp5/3.;\r
-  return rval;\r
-}\r

diff --git a/src/build_A.mexa64 b/src/build_A.mexa64
index 89a686873c4c08e9b820971886b0aede28b6eec9..ee6b1c6a38f6067822ca43508fe35aba8cbbd3fa 100755 (executable)
Binary files a/src/build_A.mexa64 and b/src/build_A.mexa64 differ

diff --git a/src/build_A2.m b/src/build_A2.m
index f6765a92be776845636ed9368776106f887f857a..9caa83afc2f1eeda6f09f6ea04b3c259103f1f67 100644 (file)
--- a/src/build_A2.m
+++ b/src/build_A2.m
@@ -4,11 +4,11 @@ function A = build_A2(coordinates,elements)
     A1 = zeros(N);
 
     % untere schranke s t obere schranke k l
-    intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ...
-        f(k1,k2,l1,l2)-f(k1,k2,l1,t2)-f(k1,k2,t1,l2)+f(k1,k2,t1,t2)...
-        -f(k1,s2,l1,l2)+f(k1,s2,l1,t2)+f(k1,s2,t1,l2)-f(k1,s2,t1,t2)...
-        -f(s1,k2,l1,l2)+f(s1,k2,l1,t2)+f(s1,k2,t1,l2)-f(s1,l2,t1,t2)...
-        +f(s1,s2,l1,l2)-f(s1,0,l1,t2)-f(s1,s2,t1,l2)+f(s1,s2,t1,t2);
+    intF = @(f,a,b,c,d,r,t,u,v)...
+        f(b,d,t,v)-f(b,d,t,u)-f(b,d,r,v)+f(b,d,r,u)...
+        -f(b,c,t,v)+f(b,c,t,u)+f(b,c,r,v)-f(b,c,r,u)...
+        -f(a,d,t,v)+f(a,d,t,u)+f(a,d,r,v)-f(a,d,r,u)...
+        +f(a,c,t,v)-f(a,c,t,u)-f(a,c,r,v)+f(a,c,r,u);
     tic
     for j = 1:N
         for k = 1:N
@@ -17,7 +17,7 @@ function A = build_A2(coordinates,elements)
             d = ek(1,:) - ej(1,:);
 %             d = ones(1,3);
             A1(j,k) = intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
-                ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));
+                ,ej(1,1),ej(2,1),ej(1,2), ej(3,2),ek(1,1), ek(2,1),ek(1,2), ek(3,2));
         end
     end
     

diff --git a/src/mex_Fpar.cpp b/src/mex_Fpar.cpp
new file mode 100644 (file)
index 0000000..9af5fac
--- /dev/null
+++ b/src/mex_Fpar.cpp
@@ -0,0 +1,44 @@
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+#include "SLPrecangle.hpp"
+
+#define my_PI 3.141592653589793
+#define EPS 0.02
+
+using namespace std;
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+       //sicherheitsabfragen zu Datengroessen
+       if (nrhs != 7)
+               mexErrMsgTxt("expected (1,x2,y1,y2,d1,d2,d3)");
+       if (nlhs > 2)
+               mexErrMsgTxt("has maximal two output arguments");
+
+
+       plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
+       plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
+
+//     plhs[0] = mxCreateDoubleMatrix(1, 2, mxREAL);
+
+
+       //printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+
+       //sol = G00(p,y1,y2,x1,x2,l);
+
+       *mxGetPr(plhs[0]) = F_par(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+       *mxGetPr(plhs[1]) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+
+//     printf("%f",FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6])));
+
+//     *(mxGetPr(plhs[0])+1) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+
+
+}
+
+
+

diff --git a/src/mex_G00.cpp b/src/mex_G00.cpp
new file mode 100644 (file)
index 0000000..0bc77f3
--- /dev/null
+++ b/src/mex_G00.cpp
@@ -0,0 +1,38 @@
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+#include "SLPrecangle.hpp"
+
+#define my_PI 3.141592653589793
+#define EPS 0.02
+
+using namespace std;
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+       //sicherheitsabfragen zu Datengroessen
+       if (nrhs != 6)
+               mexErrMsgTxt("expected (p, y1, y2, x1, x2, l)");
+       if (nlhs > 2)
+               mexErrMsgTxt("has maximal two output arguments");
+
+
+       plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
+       plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
+
+
+       //printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+
+       //sol = G00(p,y1,y2,x1,x2,l);
+
+       *mxGetPr(plhs[0]) = G00(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
+       *mxGetPr(plhs[1]) = slpADLO(*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
+
+
+}
+
+
+

diff --git a/src/mex_g0.cpp b/src/mex_g0.cpp
index b31afe499f14c1d35e7efb51dd0e9d062b5bf91d..a93de87a7be7bac614f0a4f4a784fa76e0d4e5ad 100644 (file)
--- a/src/mex_g0.cpp
+++ b/src/mex_g0.cpp
@@ -17,7 +17,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
        if (nrhs != 4)
                mexErrMsgTxt("expected (p, y, x, l)");
        if (nlhs > 2)
-               mexErrMsgTxt("has maximal two output argument");
+               mexErrMsgTxt("has maximal two output arguments");
 
 
        plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);

diff --git a/src/refineQuad.m b/src/refineQuad.m
index 8d31fc9d25829baceaa79a2dab019bb3b2e76e90..bef2205bcfa3c9d4f36338d24b36c57cb7b7002d 100644 (file)
--- a/src/refineQuad.m
+++ b/src/refineQuad.m
@@ -1,5 +1,9 @@
 function [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type)
 
+if([1 1] == size(type))
+    type = repmat(type, size(elements,1),1);
+end
+
 c_loop = size(elements,1);
 
 fa2so = repmat([1:c_loop]',1,4);

diff --git a/src/surfDoubleQuad.m b/src/surfDoubleQuad.m
new file mode 100644 (file)
index 0000000..49a9e3e
--- /dev/null
+++ b/src/surfDoubleQuad.m
@@ -0,0 +1,45 @@
+function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
+%Vierfachintegral ueber f(x1,x2,y1,y2) auf [a b]x[c d]x[r t]x[u v] mit 
+%der Simpsonregel und w Auswertungspunkten
+    s = (b-a)/w;
+    in = 0;
+    for i = a:s:b-s
+        in = in + fx2(f,i,c,d,r,t,u,v,w) + ...
+                4*fx2(f,(2*i+s)/2,c,d,r,t,u,v,w) + ...
+                  fx2(f,i+s,c,d,r,t,u,v,w);
+    end
+    in = in *s/6;
+end
+
+function in = fx2(f,x1,c,d,r,t,u,v,w)
+    s = (d-c)/w;
+    in = 0;
+    for i = c:s:d-s
+        in = in + fy1(f,x1,i,r,t,u,v,w) + ...
+                4*fy1(f,x1,(2*i+s)/2,r,t,u,v,w) + ...
+                  fy1(f,x1,i+s,r,t,u,v,w);
+    end
+    in = in *s/6;
+end
+
+function in = fy1(f,x1,x2,r,t,u,v,w)
+    s = (t-r)/w;
+    in = 0;
+    for i = r:s:t-s
+        in = in + fy2(f,x1,x2,i,u,v,w) + ...
+                4*fy2(f,x1,x2,(2*i+s)/2,u,v,w) + ...
+                  fy2(f,x1,x2,i+s,u,v,w);
+    end
+    in = in *s/6;
+end
+
+function in = fy2(f,x1,x2,y1,u,v,w)
+    s = (v-u)/w;
+    in = 0;
+    for i = u:s:v-s
+        in = in + f(x1,x2,y1,i) + ...
+                4*f(x1,x2,y1,(2*i+s)/2) + ...
+                  f(x1,x2,y1,i+s);
+    end
+    in = in *s/6;
+end
\ No newline at end of file

diff --git a/src/surfQuad.m b/src/surfQuad.m
new file mode 100644 (file)
index 0000000..b089675
--- /dev/null
+++ b/src/surfQuad.m
@@ -0,0 +1,19 @@
+function in = surfQuad(f,a,b,c,d,v)\r
+%Doppelintegral ueber f(x,y) auf [a b]x[c d] mit der Simpsonregel und v\r
+%Auswertungspunkten\r
+    s = (b-a)/v;\r
+    in = 0;\r
+    for i = a:s:b-s\r
+        in = in + fy(f,i,c,d,v) + 4*fy(f,(2*i+s)/2,c,d,v) + fy(f,i+s,c,d,v);\r
+    end\r
+    in = in *s/6;\r
+end\r
+\r
+function in = fy(f,x,c,d,v)\r
+    s = (d-c)/v;\r
+    in = 0;\r
+    for i = c:s:d-s\r
+        in = in + f(x,i) + 4*f(x,(2*i+s)/2) + f(x,i+s);\r
+    end\r
+    in = in *s/6;\r
+end
\ No newline at end of file

diff --git a/src/test_Fpar.m b/src/test_Fpar.m
new file mode 100644 (file)
index 0000000..27c3e6c
--- /dev/null
+++ b/src/test_Fpar.m
@@ -0,0 +1,19 @@
+function [ret] = test_Fpar(a,b,c,d,r,t,u,v,l1,l2,l3)
+
+
+    quad_sol = surfDoubleQuad(@(x1,x2,y1,y2)1/sqrt((x1-y1-l1).^2+(x2-y2-l2).^2+l3.^2)...
+                ,a,b,c,d,r,t,u,v,4);
+
+    intF = @(f,a,b,c,d,r,t,u,v)...
+        f(b,d,t,v)-f(b,d,t,u)-f(b,d,r,v)+f(b,d,r,u)...
+        -f(b,c,t,v)+f(b,c,t,u)+f(b,c,r,v)-f(b,c,r,u)...
+        -f(a,d,t,v)+f(a,d,t,u)+f(a,d,r,v)-f(a,d,r,u)...
+        +f(a,c,t,v)-f(a,c,t,u)-f(a,c,r,v)+f(a,c,r,u);
+            
+    sol = intF(@(x1,x2,y1,y2)mex_Fpar(x1,x2,y1,y2,l1,l2,l3),a,b,c,d,r,t,u,v);
+
+
+
+    ret = [quad_sol sol];
+end
+

diff --git a/src/test_G00.m b/src/test_G00.m
new file mode 100644 (file)
index 0000000..c23354f
--- /dev/null
+++ b/src/test_G00.m
@@ -0,0 +1,23 @@
+function [ret] = test_G00(p, x1, x2,l,a,b,c,d)
+
+
+    quad_sol = surfQuad(@(y1,y2)((x1-y1).^2+(x2-y2).^2+l.^2).^p,a,b,c,d,4);
+
+
+    [my(1) mk(1)] = mex_G00(p,a,c,x1,x2,l);
+    [my(2) mk(2)] = mex_G00(p,b,d,x1,x2,l);
+    [my(3) mk(3)] = mex_G00(p,b,c,x1,x2,l);
+    [my(4) mk(4)] = mex_G00(p,a,d,x1,x2,l);
+
+
+    my_sol = my(1)+my(2)-my(3)-my(4);
+    if(p==-0.5)
+        mk_sol = mk(1)+mk(2)-mk(3)-mk(4);
+    else
+        mk_sol = nan;
+    end
+
+
+    ret = [quad_sol my_sol mk_sol];
+end
+

diff --git a/src/test_functions.m b/src/test_functions.m
new file mode 100644 (file)
index 0000000..00ba036
--- /dev/null
+++ b/src/test_functions.m
@@ -0,0 +1,60 @@
+function done=test_functions(eps)
+
+    done = 0;
+%% g0 testen  
+p_g0 = [0.5 0 -0.5 -1 -1.5];
+
+    for p = p_g0
+        s = test_g0(p,1.4,2,1,4);
+        if(abs(s(1)-s(2))>eps)
+            p
+            s
+            return;
+        end
+    end
+    
+    % x muss ausserhalb des Intervalls liegen
+    for p = p_g0
+        s = test_g0(p,1.4,0,2,4);
+        if(abs(s(1)-s(2))>eps)
+            p
+            s
+            return;
+        end
+    end
+
+%% G00 testen
+p_G00 = [-0.5 -1.5];
+
+    for p = p_G00
+        s = test_G00(p,1.4,0.7,1,2,4,1,3);
+        if(abs(s(1)-s(2))>eps)
+            p
+            s
+            return;
+        end
+    end
+    
+    for p = p_G00
+        s = test_G00(p,1.1,0.3,0,2,4.5,1,2.5);
+        if(abs(s(1)-s(2))>eps)
+            p
+            s
+            return;
+        end
+    end
+    
+%% F_par testen
+    
+        
+        s = test_Fpar(1,3,2,3,4,5,1,3,1,3,1);
+        if(abs(s(1)-s(2))>eps)
+            p
+            s
+            return;
+        end
+        
+        
+        
+    done = 1;
+end
\ No newline at end of file

diff --git a/src/test_g0.m b/src/test_g0.m
index b08aaf032480982c18549fee5f79a64e8b49954e..1bd8d15005f16e59d3a9f75922d03e29857e3997 100644 (file)
--- a/src/test_g0.m
+++ b/src/test_g0.m
@@ -3,6 +3,10 @@ function [ret] = test_g0(p, x,l,a,b)
 
 quad_sol = quad(@(y)((x-y).^2+l.^2).^p,a,b);
 
+X = 0:0.2:5;
+Y = ((x-X).^2+l.^2).^p;
+
+plot(X,Y)
 
 [my(1) mk(1)] = mex_g0(p,a,x,l);
 [my(2) mk(2)] = mex_g0(p,b,x,l);

diff --git a/src/test_solveError.m b/src/test_solveError.m
index 7424977efacf073f996f465dcdea256490f6dae8..9a1c633e73ee018446c23f2a92f664bca1cfb5f1 100644 (file)
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -2,29 +2,36 @@
 \r
 load exmpl_2DLShape\r
 \r
+\r
+[coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
+[coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
+\r
+\r
 A = build_A(coordinates,elements);\r
+sum(sum (A-A'))\r
 b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));\r
 x = A\b\r
 xe = x'*A*x;\r
 \r
 [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
 A_fine = build_A(coordinates_fine,elements_fine);\r
+sum(sum(A_fine-A_fine'))\r
 b = sqrt(sum(quad_norm(coordinates_fine,elements_fine,'w').^2,2));\r
 x_fine = A_fine\b;\r
 xe_fine = x_fine'*A_fine*x_fine;\r
 \r
 x_interpol = zeros(size(elements_fine,1),1);\r
 for i = 1:size(elements,1)\r
-    x_interpol(f2s(i,:)) = x(i)/4;\r
+    x_interpol(f2s(i,:)) = x(i);\r
 end\r
 xe_interpol = x_interpol'*A_fine*x_interpol;\r
 \r
 [x_fine x_interpol x_fine-x_interpol]\r
 \r
-l = x_fine(f2s(1,:));\r
-T4 = [1 1 1 1; 1 -1 1 -1; 1 1 -1 -1;1 -1 -1 1];\r
-l\r
-T4*l*1/4\r
+% l = x_fine(f2s(1,:));\r
+% T4 = [1 1 1 1; 1 -1 1 -1; 1 1 -1 -1;1 -1 -1 1];\r
+% l\r
+% T4*l*1/4\r
 \r
 [xe xe_fine xe_interpol]\r
 \r

diff --git a/src/test_solveS.m b/src/test_solveS.m
index 1de19e0ee8293db082d9eb7aa47b8c870d450194..215c0c821638afae4e9eb76d03e3268409c57a0e 100644 (file)
--- a/src/test_solveS.m
+++ b/src/test_solveS.m
@@ -14,24 +14,7 @@ intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ...
     -f(s1,k2,l1,l2)+f(s1,k2,l1,t2)+f(s1,k2,t1,l2)-f(s1,l2,t1,t2)...
     +f(s1,s2,l1,l2)-f(s1,0,l1,t2)-f(s1,s2,t1,l2)+f(s1,s2,t1,t2);
 tic
-for j = 1:N
-    for k = 1:N
-        ej = coordinates(elements(j,:)',:);
-        ek = coordinates(elements(k,:)',:);
-        d = ek(1,:) - ej(1,:);
-%         fprintf('%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n'...
-%             ,j,k,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2)...
-%             ,d(1),d(2),d(3));
-        tmp = intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
-            ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));
-
-%        fprintf('%d%d|%.2f\n',j,k,tmp);
-        A1(j,k) = 1/ (4*pi) *tmp;
-%           intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
-%             ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));
-
-    end
-end
+A2 = build_A2(coordinates,elements);
 t(1) = toc;
 disp ' '
 tic