From ccd1d4c4b3d79205eed91f67c0ee3fdc9e4e39f3 Mon Sep 17 00:00:00 2001
From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Fri, 20 May 2011 21:45:17 +0000
Subject: [PATCH] =?utf8?q?[src]=20Fehler=20gefunde:=20Lag=20in=20den=20Gre?=
 =?utf8?q?nzen=20[src]=20neue=20=C3=9Cberpr=C3=BCfung=20f=C3=BCr=20die=20L?=
 =?utf8?q?age=20von=20parallenen=20Rechtecken?=
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@27 26120e32-c555-405d-b3e1-1f783fb42516
---
 src/SLPrecangle.cpp   |  11 +++++
 src/build_A.cpp       | 108 ++++++++++++++++++++++++++++++++----------
 src/build_A2.m        |  26 ++++++----
 src/surfDoubleQuad.m  |   2 +-
 src/test_solveError.m |  21 ++++----
 5 files changed, 122 insertions(+), 46 deletions(-)

diff --git a/src/SLPrecangle.cpp b/src/SLPrecangle.cpp
index 4a67229..42b08e4 100644
--- a/src/SLPrecangle.cpp
+++ b/src/SLPrecangle.cpp
@@ -277,6 +277,17 @@ double quadInt(
 			*/
 }
 
+double quad0Int(
+		double(*f)(double, double, double, double, double, double, double),
+		double b, double d, double t, double v, double d1, double d2, double d3) {
+
+   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)
+    -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)
+    -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)
+    +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);
+
+}
+
 
 
 
diff --git a/src/build_A.cpp b/src/build_A.cpp
index 6050480..ee198db 100644
--- a/src/build_A.cpp
+++ b/src/build_A.cpp
@@ -4,6 +4,9 @@
 #include "mex.h"
 #include <stdlib.h>
 
+//#include "tbb/parallel_for.h"
+
+
 #include "SLPrecangle.hpp"
 
 #define my_PI 3.141592653589793
@@ -12,6 +15,9 @@
 using namespace std;
 
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+
+	int i,j,k;
+
 	//sicherheitsabfragen zu Datengroessen
 	if (nrhs != 2)
 		mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx3))");
@@ -38,20 +44,25 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
 	double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
 	double * xn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
 
 	double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
 	double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
 	double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
 	double * yn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
 
 	double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
 
 	double tmp;
 
-	int rx, ry;
+	//LageInformationen
+	int rx, rxa, rxb, ry, rya, ryb;
 
 	//Ausrechnen
-	for (int j = 0; j < em; ++j) {
+	for (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];
@@ -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];
 		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]);
+		for (i = 0; i<3;++i)
+			xa[i] = x2[i] - x1[i];
+
+		for (i = 0; i<3;++i)
+			xb[i] = x3[i] - x1[i];
+
+		//Kreuzprodukt
+		xn[0] = (xa[1]) * (xb[2]) - (xa[2]) * (xb[1]);
+		xn[1] = (xa[2]) * (xb[0]) - (xa[0]) * (xb[2]);
+		xn[2] = (xa[0]) * (xb[1]) - (xa[1]) * (xb[0]);
+
+
 
 		//if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])
 		if (xn[2] != 0)
@@ -79,9 +96,25 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 		else
 			rx = 0;
 
+		if (xa[2] != 0)
+			rxa = 2;
+		else if (xa[1] != 0)
+			rxa = 1;
+		else
+			rxa = 0;
+
+		if (xb[2] != 0)
+			rxb = 2;
+		else if (xb[1] != 0)
+			rxb = 1;
+		else
+			rxb = 0;
+
+
+
 		// printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);
 
-		for (int k = 0; k < em; ++k) {
+		for (k = 0; k < em; ++k) {
 			y1[0] = C[(int) E[k] - 1];
 			y1[1] = C[cm + (int) E[k] - 1];
 			y1[2] = C[2 * cm + (int) E[k] - 1];
@@ -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];
 			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]);
+			for (i = 0; i<3;++i)
+				ya[i] = y2[i] - y1[i];
+
+			for (i = 0; i<3;++i)
+				yb[i] = y3[i] - y1[i];
+
+			yn[0] = (ya[1]) * (yb[2]) - (ya[2]) * (yb[1]);
+			yn[1] = (ya[2]) * (yb[0]) - (ya[0]) * (yb[2]);
+			yn[2] = (ya[0]) * (yb[1]) - (ya[1]) * (yb[0]);
 
 			if (yn[2] != 0)
 				ry = 2;
@@ -108,24 +144,44 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 			else
 				ry = 0;
 
-			d[0] = y1[0] - x1[0];
-			d[1] = y1[1] - x1[1];
-			d[2] = y1[2] - x1[2];
+			if (xa[2] != 0)
+				rxa = 2;
+			else if (xa[1] != 0)
+				rxa = 1;
+			else
+				rxa = 0;
+
+			if (xb[2] != 0)
+				rxb = 2;
+			else if (xb[1] != 0)
+				rxb = 1;
+			else
+				rxb = 0;
+
+			for (i = 0; i<3;++i)
+				d[i] = y1[i] - x1[i];
 
-			if (rx = ry) {
-				if (rx = 2) {
+			printf("(%d,%d)",rx,ry);
+			if (rx == ry) {
+				if (rxa == rya) {
+					//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],
+					//		y2[0], y1[1], y3[1], d[0], d[1], d[2]);
+					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
+					tmp = quadInt(F_par, 0, xa[rxa], 0, xb[rxb], 0,
+							ya[rxa], 0, yb[rxb], d[rxa], d[rxb], d[rx]);
+					//	 printf("%d%d|%.2f\n",j,k,tmp);
+				} else{
 					//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],
 					//		y2[0], y1[1], y3[1], d[0], d[1], d[2]);
 					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
-					tmp = quadInt(F_par, x1[0], x2[0], x1[1], x3[1], y1[0],
-							y2[0], y1[1], y3[1], d[0], d[1], d[2]);
+					tmp = quadInt(F_par, 0, xa[rxa], 0, xb[rxb], 0,
+							yb[rxa], 0, ya[rxb], d[rxa], d[rxb], d[rx]);
 					//	 printf("%d%d|%.2f\n",j,k,tmp);
-					A[(k * em) + j] = 1. / (4 * my_PI) * tmp;
+				}
+				A[(k * em) + j] = 1. / (4 * my_PI) * tmp;
 
-				} else
-					A[(k * em) + j] = 0;
 			} else {
-				A[(k * em) + j] = 0;
+				A[(k * em) + j] = NAN;
 			}
 
 			// Vorbereiten der DATEN
diff --git a/src/build_A2.m b/src/build_A2.m
index 4956416..9550b53 100644
--- 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 757d135..a93accd 100644
--- 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 366ea25..3a9fba0 100644
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -1,23 +1,26 @@
 
 
 load exmpl_2DLShape
-%[coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
+
+coordinates = coordinates(:,[ 1 3 2]);
+
+% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
+% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
- [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
 % build_A(coordinates,elements);
 
-A = build_A2(coordinates,elements);
-sum(sum (A-A'));
-b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2))
-x = A\b;
-xe = x'*A*x
+% A = build_A2(coordinates,elements);
+% sum(sum (A-A'));
+% b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));
+% x = A\b;
+% xe = x'*A*x
 
-A = build_A(coordinates,elements);% .* (ones(size(elements,1))+eye(size(elements,1)));
+A = build_A(coordinates,elements)
 sum(sum (A-A'));
 b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));
-x = A\b
+x = A\b;
 xe = x'*A*x
 
 % [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
-- 
2.47.3