[src] mex File läuft durch und sollte das Richtige machen

[src] Fehlerhafte Funktionen ! (Übergebene Parameter stimmen, kommen aber auf verschiedene Lösungen)

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

diff --git a/src/build_A.cpp b/src/build_A.cpp
index 021f5b65be0eab85276af5e405bb44e76160eb9d..3c4e9fbb48ada9c8d0f6fff39e708d838acbee57 100644 (file)
--- a/src/build_A.cpp
+++ b/src/build_A.cpp
@@ -65,22 +65,93 @@ void mexFunction(
     double * yn = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL));\r
     \r
     double * d = mxGetPr(mxCreateDoubleMatrix(3,1,mxREAL));\r
+\r
+       double tmp;\r
+\r
+\r
+       int rx,ry;\r
     \r
     \r
     \r
     //Ausrechnen\r
 for (int 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
+\r
+        x2[0] = C[(int)E[em+j]-1];\r
+        x2[1] = C[cm+(int)E[em+j]-1];\r
+        x2[2] = C[2*cm+(int)E[em+j]-1];\r
+\r
+        x3[0] = C[(int)E[2*em+j]-1];\r
+        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]-x1[1]);\r
+        xn[1] = (x2[2]-x1[2])*(x3[0]-x1[0])-(x2[0]-x1[0])*(x3[2]-x1[2]);\r
+        xn[2] = (x2[0]-x1[0])*(x3[1]-x1[1])-(x2[1]-x1[1])*(x3[0]-x1[0]);\r
+\r
+        //if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])\r
+        if(xn[2]!=0)\r
+                rx = 2;\r
+        else if(xn[1]!=0)\r
+                rx = 1;\r
+        else\r
+                rx = 0;\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
-        \r
-     x1[0] = E(j);\r
-        \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
+\r
+        y2[0] = C[(int)E[em+k]-1];\r
+        y2[1] = C[cm+(int)E[em+k]-1];\r
+        y2[2] = C[2*cm+(int)E[em+k]-1];\r
+\r
+        y3[0] = C[(int)E[2*em+k]-1];\r
+        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])*(y3[1]-y1[1]);\r
+        yn[1] = (y2[2]-y1[2])*(y3[0]-y1[0])-(y2[0]-y1[0])*(y3[2]-y1[2]);\r
+        yn[2] = (y2[0]-y1[0])*(y3[1]-y1[1])-(y2[1]-y1[1])*(y3[0]-y1[0]);\r
+\r
+        if(yn[2]!=0)\r
+                ry = 2;\r
+        else if(yn[1]!=0)\r
+                ry = 1;\r
+        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
+\r
+        if(rx=ry){\r
+                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
+                        tmp = quadInt(F_par,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|%.2f\n",j,k,tmp);\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
+        }\r
+\r
+\r
+\r
+ \r
      // Vorbereiten der DATEN\r
      //   ej = coordinates(elements(j,:)',:);\r
      //   ek = coordinates(elements(k,:)',:);\r
      //   d = ek(1,:) - ej(1,:);\r
 \r
      // AUSRECHNEN\r
-        A[(k*em)+j] = 1./(4*my_PI);// *\r
+  //     A[(k*em)+j] = 1./(4*my_PI);// *\r
          // quadInt(F_par,\r
          //   ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));\r
 \r
@@ -89,7 +160,15 @@ for (int j=0;j<em;++j){
    \r
 }    \r
      //Rueckgabe (eventuell zurueck schreiben)     \r
-    \r
+ //   mxFree(x1);\r
+       //mxFree(x2);\r
+       //mxFree(x3);\r
+       //mxFree(y1);\r
+       //mxFree(y2);\r
+       //mxFree(y3);\r
+       //mxFree(xn);\r
+       //mxFree(yn);\r
+       //mxFree(d);\r
     \r
  return;   \r
 }\r

diff --git a/src/test_solve.m b/src/test_solve.m
index 99c5f2327dc2d18e4a4f028af5f83ae7cf49f5a8..b40f953d0e1d9f711cfd994e7562906be75a61cc 100644 (file)
--- a/src/test_solve.m
+++ b/src/test_solve.m
@@ -1,7 +1,10 @@
 
-N = size(elements);
+N = size(elements,1);
 
-A = zeros(N);
+A1 = zeros(N);
+
+x = zeros(N,2);
+xe = zeros(1,2);
 
 % untere schranke s t obere schranke k l
 intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ...
@@ -15,14 +18,30 @@ for j = 1:N
         ej = coordinates(elements(j,:)',:);
         ek = coordinates(elements(k,:)',:);
         d = ek(1,:) - ej(1,:);
-
-        A(j,k) = 1/ (4*pi) *...
-          intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
+%         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
+disp ' '
+
+A2 = build_A(coordinates,elements);
+
+b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));
+
+x(:,1) = A1\b;
+xe(1) = x(:,1)'*A1*x(:,1);
+
+x(:,2) = A2\b;
+xe(2) = x(:,2)'*A2*x(:,2);
 
-A
-b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2))
-x = A\b
\ No newline at end of file
+x
+xe
\ No newline at end of file