From 696e89c1007796e37c75bc6b99091d08fc82974e Mon Sep 17 00:00:00 2001
From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Thu, 21 Apr 2011 15:45:00 +0000
Subject: [PATCH] =?utf8?q?[src]=20mex=20File=20l=C3=A4uft=20durch=20und=20?=
 =?utf8?q?sollte=20das=20Richtige=20machen=20[src]=20Fehlerhafte=20Funktio?=
 =?utf8?q?nen=20!=20(=C3=9Cbergebene=20Parameter=20stimmen,=20kommen=20abe?=
 =?utf8?q?r=20auf=20verschiedene=20L=C3=B6sungen)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=utf8
Content-Transfer-Encoding: 8bit

git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@12 26120e32-c555-405d-b3e1-1f783fb42516
---
 src/build_A.cpp  | 89 +++++++++++++++++++++++++++++++++++++++++++++---
 src/test_solve.m | 35 ++++++++++++++-----
 2 files changed, 111 insertions(+), 13 deletions(-)

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