[src] Fehler gefunde: Lag in den Grenzen

[src] neue Überprüfung für die Lage von parallenen Rechtecken

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

diff --git a/src/SLPrecangle.cpp b/src/SLPrecangle.cpp
index 4a6722947d58b211a6cb834ff4abee78b540ed6e..42b08e46b661e11c2de65cdabc1f60705307330e 100644 (file)
--- a/src/SLPrecangle.cpp
+++ b/src/SLPrecangle.cpp
@@ -277,6 +277,17 @@ double quadInt(
                        */\r
 }\r
 \r
+double quad0Int(\r
+               double(*f)(double, double, double, double, double, double, double),\r
+               double b, double d, double t, double v, double d1, double d2, double d3) {\r
+\r
+   return f(b,d,t,v,d1,d2,d3)-f(b,d,t,0,d1,d2,d3)-f(b,d,0,v,d1,d2,d3)+f(b,d,0,0,d1,d2,d3)\r
+    -f(b,0,t,v,d1,d2,d3)+f(b,0,t,0,d1,d2,d3)+f(b,0,0,v,d1,d2,d3)-f(b,0,0,0,d1,d2,d3)\r
+    -f(0,d,t,v,d1,d2,d3)+f(0,d,t,0,d1,d2,d3)+f(0,d,0,v,d1,d2,d3)-f(0,d,0,0,d1,d2,d3)\r
+    +f(0,0,t,v,d1,d2,d3)-f(0,0,t,0,d1,d2,d3)-f(0,0,0,v,d1,d2,d3)+f(0,0,0,0,d1,d2,d3);\r
+\r
+}\r
+\r
 \r
 \r
 \r

diff --git a/src/build_A.cpp b/src/build_A.cpp
index 60504802cdeda10b9594b8a88db745a3f436c16b..ee198dbdd426f2d23b164c143b8f8a1a7f89492c 100644 (file)
--- a/src/build_A.cpp
+++ b/src/build_A.cpp
@@ -4,6 +4,9 @@
 #include "mex.h"\r
 #include <stdlib.h>\r
 \r
+//#include "tbb/parallel_for.h"\r
+\r
+\r
 #include "SLPrecangle.hpp"\r
 \r
 #define my_PI 3.141592653589793\r
@@ -12,6 +15,9 @@
 using namespace std;\r
 \r
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
+\r
+       int i,j,k;\r
+\r
        //sicherheitsabfragen zu Datengroessen\r
        if (nrhs != 2)\r
                mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx3))");\r
@@ -38,20 +44,25 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
        double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
        double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
        double * xn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
 \r
        double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
        double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
        double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
        double * yn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
 \r
        double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
 \r
        double tmp;\r
 \r
-       int rx, ry;\r
+       //LageInformationen\r
+       int rx, rxa, rxb, ry, rya, ryb;\r
 \r
        //Ausrechnen\r
-       for (int j = 0; j < em; ++j) {\r
+       for (j = 0; j < em; ++j) {\r
                x1[0] = C[(int) E[j] - 1];\r
                x1[1] = C[cm + (int) E[j] - 1];\r
                x1[2] = C[2 * cm + (int) E[j] - 1];\r
@@ -64,12 +75,18 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                x3[1] = C[cm + (int) E[2 * em + j] - 1];\r
                x3[2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
 \r
-               xn[0] = (x2[1] - x1[1]) * (x3[2] - x1[2]) - (x2[2] - x1[2]) * (x3[1]\r
-                               - x1[1]);\r
-               xn[1] = (x2[2] - x1[2]) * (x3[0] - x1[0]) - (x2[0] - x1[0]) * (x3[2]\r
-                               - x1[2]);\r
-               xn[2] = (x2[0] - x1[0]) * (x3[1] - x1[1]) - (x2[1] - x1[1]) * (x3[0]\r
-                               - x1[0]);\r
+               for (i = 0; i<3;++i)\r
+                       xa[i] = x2[i] - x1[i];\r
+\r
+               for (i = 0; i<3;++i)\r
+                       xb[i] = x3[i] - x1[i];\r
+\r
+               //Kreuzprodukt\r
+               xn[0] = (xa[1]) * (xb[2]) - (xa[2]) * (xb[1]);\r
+               xn[1] = (xa[2]) * (xb[0]) - (xa[0]) * (xb[2]);\r
+               xn[2] = (xa[0]) * (xb[1]) - (xa[1]) * (xb[0]);\r
+\r
+\r
 \r
                //if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])\r
                if (xn[2] != 0)\r
@@ -79,9 +96,25 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                else\r
                        rx = 0;\r
 \r
+               if (xa[2] != 0)\r
+                       rxa = 2;\r
+               else if (xa[1] != 0)\r
+                       rxa = 1;\r
+               else\r
+                       rxa = 0;\r
+\r
+               if (xb[2] != 0)\r
+                       rxb = 2;\r
+               else if (xb[1] != 0)\r
+                       rxb = 1;\r
+               else\r
+                       rxb = 0;\r
+\r
+\r
+\r
                // printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);\r
 \r
-               for (int k = 0; k < em; ++k) {\r
+               for (k = 0; k < em; ++k) {\r
                        y1[0] = C[(int) E[k] - 1];\r
                        y1[1] = C[cm + (int) E[k] - 1];\r
                        y1[2] = C[2 * cm + (int) E[k] - 1];\r
@@ -94,12 +127,15 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                        y3[1] = C[cm + (int) E[2 * em + k] - 1];\r
                        y3[2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
 \r
-                       yn[0] = (y2[1] - y1[1]) * (y3[2] - y1[2]) - (y2[2] - y1[2])\r
-                                       * (y3[1] - y1[1]);\r
-                       yn[1] = (y2[2] - y1[2]) * (y3[0] - y1[0]) - (y2[0] - y1[0])\r
-                                       * (y3[2] - y1[2]);\r
-                       yn[2] = (y2[0] - y1[0]) * (y3[1] - y1[1]) - (y2[1] - y1[1])\r
-                                       * (y3[0] - y1[0]);\r
+                       for (i = 0; i<3;++i)\r
+                               ya[i] = y2[i] - y1[i];\r
+\r
+                       for (i = 0; i<3;++i)\r
+                               yb[i] = y3[i] - y1[i];\r
+\r
+                       yn[0] = (ya[1]) * (yb[2]) - (ya[2]) * (yb[1]);\r
+                       yn[1] = (ya[2]) * (yb[0]) - (ya[0]) * (yb[2]);\r
+                       yn[2] = (ya[0]) * (yb[1]) - (ya[1]) * (yb[0]);\r
 \r
                        if (yn[2] != 0)\r
                                ry = 2;\r
@@ -108,24 +144,44 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                        else\r
                                ry = 0;\r
 \r
-                       d[0] = y1[0] - x1[0];\r
-                       d[1] = y1[1] - x1[1];\r
-                       d[2] = y1[2] - x1[2];\r
+                       if (xa[2] != 0)\r
+                               rxa = 2;\r
+                       else if (xa[1] != 0)\r
+                               rxa = 1;\r
+                       else\r
+                               rxa = 0;\r
+\r
+                       if (xb[2] != 0)\r
+                               rxb = 2;\r
+                       else if (xb[1] != 0)\r
+                               rxb = 1;\r
+                       else\r
+                               rxb = 0;\r
+\r
+                       for (i = 0; i<3;++i)\r
+                               d[i] = y1[i] - x1[i];\r
 \r
-                       if (rx = ry) {\r
-                               if (rx = 2) {\r
+                       printf("(%d,%d)",rx,ry);\r
+                       if (rx == ry) {\r
+                               if (rxa == rya) {\r
+                                       //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,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 %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
+                                       tmp = quadInt(F_par, 0, xa[rxa], 0, xb[rxb], 0,\r
+                                                       ya[rxa], 0, yb[rxb], d[rxa], d[rxb], d[rx]);\r
+                                       //       printf("%d%d|%.2f\n",j,k,tmp);\r
+                               } else{\r
                                        //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,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 %.1f %.1f %.1f\n",j,k,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
+                                       tmp = quadInt(F_par, 0, xa[rxa], 0, xb[rxb], 0,\r
+                                                       yb[rxa], 0, ya[rxb], d[rxa], d[rxb], d[rx]);\r
                                        //       printf("%d%d|%.2f\n",j,k,tmp);\r
-                                       A[(k * em) + j] = 1. / (4 * my_PI) * tmp;\r
+                               }\r
+                               A[(k * em) + j] = 1. / (4 * my_PI) * tmp;\r
 \r
-                               } else\r
-                                       A[(k * em) + j] = 0;\r
                        } else {\r
-                               A[(k * em) + j] = 0;\r
+                               A[(k * em) + j] = NAN;\r
                        }\r
 \r
                        // Vorbereiten der DATEN\r

diff --git a/src/build_A2.m b/src/build_A2.m
index 4956416d102d37a5e2109d8d2389cda33168913c..9550b53dc8464b93f65b7ce3445dee8d7967da97 100644 (file)
--- a/src/build_A2.m
+++ b/src/build_A2.m
@@ -4,23 +4,29 @@ function A = build_A2(coordinates,elements)
     A1 = zeros(N);
 
     % untere schranke s t obere schranke k l
-%     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
+    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);
+    
     for j = 1:N
         for k = 1:N
             ej = coordinates(elements(j,:)',:);
             ek = coordinates(elements(k,:)',:);
+            
             d = ek(1,:) - ej(1,:);
+            
+            ej = ej - repmat(ej(1,:),3,1);
+            ek = ek - repmat(ek(1,:),3,1);
+
 %             d = ones(1,3);
-            A1(j,k) = surfDoubleQuad(@(x1,x2,y1,y2) 1/sqrt((x1-y1-d(1)).^2+(x2-y2-d(2)).^2+d(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),4);
 
-%             A1(j,k) = intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
-%                 ,ej(1,1),ej(2,1),ej(1,2), ej(3,2),ek(1,1), ek(2,1),ek(1,2), ek(3,2));
+%             A1(j,k) = surfDoubleQuad(@(x1,x2,y1,y2) 1/sqrt((x1-y1-d(1)).^2+(x2-y2-d(2)).^2+d(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),4);
+
+            A1(j,k) = intF(@(x1,x2,y1,y2) mex_Fpar(x1,x2,y1,y2,d(1),d(2),d(3))...
+                ,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/surfDoubleQuad.m b/src/surfDoubleQuad.m
index 757d135497a9dccf009e092b136f651c3225bccb..a93accdc8ac85fbeb7b6b210cdf7f1759795d839 100644 (file)
--- a/src/surfDoubleQuad.m
+++ b/src/surfDoubleQuad.m
@@ -99,6 +99,6 @@ function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
                     
                 end
             end
-            in
+%             in
 
 end
\ No newline at end of file

diff --git a/src/test_solveError.m b/src/test_solveError.m
index 366ea259b910d9e389579c4a8a2c5b953759342b..3a9fba01840acc56d7e665aab48d468fcc579876 100644 (file)
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -1,23 +1,26 @@
 \r
 \r
 load exmpl_2DLShape\r
-%[coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
+\r
+coordinates = coordinates(:,[ 1 3 2]);\r
+\r
+% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
+% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
- [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
 % build_A(coordinates,elements);\r
 \r
-A = build_A2(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
+% A = build_A2(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
-A = build_A(coordinates,elements);% .* (ones(size(elements,1))+eye(size(elements,1)));\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
+x = A\b;\r
 xe = x'*A*x\r
 \r
 % [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r