From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Sun, 5 Jun 2011 19:55:24 +0000 (+0000)
Subject: [src] F_ort implementiert
X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=a82f2dd13fcfdcf41343e5c72e6321ab21a1c6ca;p=bacc.git

[src] F_ort implementiert
[src] test dazu geschrieben
[src] build_A für orthogonale Stücke erweitert
[src] surfDoubleQuad (suboptimal)
[src] example 3dFichCube gebaut [-1 1]^3
[src] code jetzt auch in Mayer eingebaut... (läuft noch nicht)
[dir] recangle -> rectangle in MayerCode

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

diff --git a/src/SLPrecangle.cpp b/src/SLPrecangle.cpp
index 42b08e4..a5f3b7e 100644
--- a/src/SLPrecangle.cpp
+++ b/src/SLPrecangle.cpp
@@ -14,158 +14,14 @@ using namespace std;
 int sign(double);
 //double arsinh(double);
 
-/*
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
-	//sicherheitsabfragen zu Datengroessen
-	if (nrhs != 2)
-		mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx3))");
-	if (nlhs > 1)
-		mexErrMsgTxt("has only one output argument");
-
-	int cm = mxGetM(prhs[0]);
-	int cn = mxGetN(prhs[0]);
-
-	if (cn != 3)
-		mexErrMsgTxt("expected coordinates (Nx3)");
-	int em = mxGetM(prhs[1]);
-	int en = mxGetN(prhs[1]);
-	if (en != 3)
-		mexErrMsgTxt("expected elements (Mx3)");
-	//Vorbereitung der Daten
-
-	plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);
-	double * A = mxGetPr(plhs[0]);
-	double * C = mxGetPr(prhs[0]);
-	double * E = mxGetPr(prhs[1]);
-
-	double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
-	double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
-	double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
-	double * xn = 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 * 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) {
-			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]);
-					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
-					tmp = quadInt(FLO_plane, 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);// *
-			// 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));
-
-		}
-
-	}
-	//Rueckgabe (eventuell zurueck schreiben)
-	//   mxFree(x1);
-	//mxFree(x2);
-	//mxFree(x3);
-	//mxFree(y1);
-	//mxFree(y2);
-	//mxFree(y3);
-	//mxFree(xn);
-	//mxFree(yn);
-	//mxFree(d);
-
-	return;
-}
-*/
 int inline sign(double x) {
 	return x > 0 ? 1 : (x < 0 ? -1 : 0);
 }
 /*
-double inline arsinh(double x) {
-	return log(x + sqrt(x * x + 1));
-}
-*/
+ double inline arsinh(double x) {
+ return log(x + sqrt(x * x + 1));
+ }
+ */
 //y-x muss != 0 sein
 double g0(double p, double y, double x, double l) {
 	//printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);
@@ -182,7 +38,7 @@ double g0(double p, double y, double x, double l) {
 		else if (p == -0.5)
 			sol = asinh((y - x) / fabs(l));
 		else if (p == -1)
-			sol = atan((y - x) / fabs(l))/l;
+			sol = atan((y - x) / fabs(l)) / l;
 		else if (p == -1.5)
 			sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l));
 		else
@@ -200,9 +56,9 @@ double g0(double p, double y, double x, double l) {
 
 double G00(double p, double y1, double y2, double x1, double x2, double l) {
 	//	printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);
-
+	double pt = p + 1.5;
 	double sol = 0;
-	if (p == -1.5) {
+	if (pt == 0) {
 		if (l == 0) {
 			sol = -sqrt((y1 - x1) * (y1 - x1) + (y2 - x2) * (y2 - x2)) / ((y1
 					- x1) * (y2 - x2));
@@ -213,7 +69,7 @@ double G00(double p, double y1, double y2, double x1, double x2, double l) {
 							+ l * l)) + 1);
 		}
 
-	} else if (p == -0.5) {
+	} else if ((pt > 0) && ((int) pt == pt)) {
 		if (l != 0)
 			sol = 2 * p * l * l * G00(p - 1, y1, y2, x1, x2, l);
 		if ((y1 - x1) != 0)
@@ -225,7 +81,7 @@ double G00(double p, double y1, double y2, double x1, double x2, double l) {
 		sol /= 2 * p + 2;
 	} else {
 		sol = NAN;
-		printf("warning in G00: no case for p defined\n");
+		printf("warning in G00: no case for p=%.2f defined\n", p);
 		//mexErrMsgTxt("no case for p defined\n");
 	}
 
@@ -239,7 +95,7 @@ double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
 	double sol = (x1 - y1 - d1) * (x2 - y2 - d2);
 
 	if (sol != 0)
-		sol *= G00(-0.5,x1, x2, y1 + d1, y2 + d2, d3);
+		sol *= G00(-0.5, x1, x2, y1 + d1, y2 + d2, d3);
 
 	if ((x1 - y1 - d1) != 0)
 		sol -= (x1 - y1 - d1) * g0(0.5, x1, y1 + d1,
@@ -255,43 +111,70 @@ double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
 	return sol;
 }
 
+double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,
+		double d3) {
+
+	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y2,y3,d1,d2,d3);
+	double sol = -G00(0.5, y3, x1, -d3, d1, x2 - y2 - d2);
+
+	if ((x1 - d1) * (x2 - y2 - d2) != 0)
+		sol -= (x1 - d1) * (x2 - y2 - d2) * G00(-0.5, x2, y3, y2 + d2, -d3,
+				x1 - d1);
+
+	if ((x1 - d1) != 0)
+		sol += (x1 - d1) * g0(0.5, y3, -d3,
+				sqrt((x1 - d1) * (x1 - d1) + (x2 - y2 - d2) * (x2 - y2 - d2)));
+
+	if ((y3 + d3) * (x2 - y2 - d2) != 0)
+		sol -= (y3 + d3) * (x2 - y2 - d2) * G00(-0.5, x1, x2, d1, y2 + d2,
+				-y3 - d3);
+
+	if ((y3 + d3) != 0)
+		sol += (y3 + d3) * g0(0.5, x1, d1,
+				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + (y3 + d3) * (y3 + d3)));
+
+	return sol / 2.;
+}
+
 double quadInt(
 		double(*f)(double, double, double, double, double, double, double),
-		double a, double b, double c, double d, double r, double t,
-		double u, double v, double d1, double d2, double d3) {
-
-   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)
-    -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)
-    -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)
-    +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);
-/*
-	return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f(
-			k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1,
-			s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2,
-			t1, l2, d1, d2, d3) - f(k1, s2, t1, t2, d1, d2, d3) - f(s1, k2, l1,
-			l2, d1, d2, d3) + f(s1, k2, l1, t2, d1, d2, d3) + f(s1, k2, t1, l2,
-			d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1,
-			d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2,
-			d3) + f(s1, s2, t1, t2, d1, d2, d3);
-
-			*/
+		double a, double b, double c, double d, double r, double t, double u,
+		double v, double d1, double d2, double d3) {
+
+	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) - 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) - 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)
+			+ 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);
+	/*
+	 return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f(
+	 k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1,
+	 s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2,
+	 t1, l2, d1, d2, d3) - f(k1, s2, t1, t2, d1, d2, d3) - f(s1, k2, l1,
+	 l2, d1, d2, d3) + f(s1, k2, l1, t2, d1, d2, d3) + f(s1, k2, t1, l2,
+	 d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1,
+	 d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2,
+	 d3) + f(s1, s2, t1, t2, d1, d2, d3);
+
+	 */
 }
 
 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);
+	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);
 
 }
 
-
-
-
-
 double slpADLO(double y1, double y2, double x1, double x2, double a) {
 	double G3 = 0;
 	double gL = 0;
@@ -354,7 +237,7 @@ double compute_g0(double p, double y, double x, double a) {
 			return (y - x) * pow((x * x + y * y + a * a - 2 * x * y), -0.5)
 					/ (a * a);
 		default:
-			printf("p must be either 0, -1/2, -1 or -3/2.");
+			printf("p must be either 0, -1/2, -1 or -3/2. (%.2f)\n", p);
 			return NAN;
 		}
 	}
@@ -375,3 +258,19 @@ double FLO_plane(double x1, double x2, double y1, double y2, double delta1,
 			x1, yd1, tmp3) - tmp2 * compute_g0(-0.5, x2, yd2, tmp4) + tmp5 / 3.;
 	return rval;
 }
+
+double FLO_perp(double x1, double x2, double y2, double y3, double delta1,
+		double delta2, double a) {
+	double xd1 = x1 - delta1;
+	double xd2 = x2 - delta2;
+	double yd2 = y2 - delta2;
+	double yd3 = y3 - a;
+	double tmp2 = x2 - y2 - delta2;
+	double rval = 0;
+
+	rval = xd1 * compute_g0(y3, a, sqrt(xd1 * xd1 + tmp2 * tmp2), -0.5) + yd3
+			* compute_g0(x1, delta1, sqrt(tmp2 * tmp2 + yd3 * yd3), -0.5) - xd1
+			* tmp2 * slpADLO(x2, y3, yd2, a, xd1) - yd3 * tmp2 * slpADLO(x1,
+			x2, delta1, yd2, -y3 - a) - slpADLO(y3, x1, -a, delta1, tmp2);
+	return rval;
+}
diff --git a/src/SLPrecangle.hpp b/src/SLPrecangle.hpp
index 418bf5e..d801cfe 100644
--- a/src/SLPrecangle.hpp
+++ b/src/SLPrecangle.hpp
@@ -10,13 +10,14 @@ double g0(double, double, double, double);
 double G00(double, double, double, double, double, double);
 // sol = F_par(x1,x2,y1,y2,d1,d2,d3);
 double F_par(double, double, double, double, double, double, double);
-// sol = F_ort(x1,x2,y1,y2,d1,d2,d3);
-//double F_ort(double, double, double, double, double, double, double);
+// sol = F_ort(x1,x2,y2,y3,d1,d2,d3);
+double F_ort(double, double, double, double, double, double, double);
 
 
 double slpADLO(double, double, double, double, double);
 double compute_g0(double, double, double, double);
 double FLO_plane(double, double, double, double, double, double, double);
+double FLO_perp(double, double, double, double, double, double, double);
 
 // sol = quadInt((F_par/F_ort),a,b,c,d,r,t,u,v,d1,d2,d3);
 double quadInt(
diff --git a/src/bem3d/src/boundary_mesh/testhilbertmesh b/src/bem3d/src/boundary_mesh/testhilbertmesh
index 51b8ad8..d49c645 100755
Binary files a/src/bem3d/src/boundary_mesh/testhilbertmesh and b/src/bem3d/src/boundary_mesh/testhilbertmesh differ
diff --git a/src/bem3d/src/laplace/SLPrecangle.cpp b/src/bem3d/src/laplace/SLPrecangle.cpp
deleted file mode 100644
index 18f174b..0000000
--- a/src/bem3d/src/laplace/SLPrecangle.cpp
+++ /dev/null
@@ -1,357 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include <stdlib.h>
-
-#include "SLPrecangle.hpp"
-
-#define my_PI 3.141592653589793
-#define EPS 0.02
-
-using namespace std;
-
-int sign(double);
-double arsinh(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);
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
-	//sicherheitsabfragen zu Datengroessen
-	if (nrhs != 2)
-		mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx3))");
-	if (nlhs > 1)
-		mexErrMsgTxt("has only one output argument");
-
-	int cm = mxGetM(prhs[0]);
-	int cn = mxGetN(prhs[0]);
-
-	if (cn != 3)
-		mexErrMsgTxt("expected coordinates (Nx3)");
-	int em = mxGetM(prhs[1]);
-	int en = mxGetN(prhs[1]);
-	if (en != 3)
-		mexErrMsgTxt("expected elements (Mx3)");
-	//Vorbereitung der Daten
-
-	plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);
-	double * A = mxGetPr(plhs[0]);
-	double * C = mxGetPr(prhs[0]);
-	double * E = mxGetPr(prhs[1]);
-
-	double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
-	double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
-	double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
-	double * xn = 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 * 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) {
-			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]);
-					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
-					tmp = quadInt(FLO_plane, 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);// *
-			// 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));
-
-		}
-
-	}
-	//Rueckgabe (eventuell zurueck schreiben)
-	//   mxFree(x1);
-	//mxFree(x2);
-	//mxFree(x3);
-	//mxFree(y1);
-	//mxFree(y2);
-	//mxFree(y3);
-	//mxFree(xn);
-	//mxFree(yn);
-	//mxFree(d);
-
-	return;
-}
-
-int inline sign(double x) {
-	return x > 0 ? 1 : (x < 0 ? -1 : 0);
-}
-
-double inline arsinh(double x) {
-	return log(x + sqrt(x * x + 1));
-}
-
-//y-x muss != 0 sein
-double g0(double p, double y, double x, double l) {
-	//printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);
-
-	double sol = 0;
-
-	if (l != 0) {
-		if (p == 0.5) {
-			sol = (y - x) / 2 * sqrt((y - x) * (y - x) + l * l) + l * l / 2
-					* arsinh((y - x) / fabs(l));
-			//          printf("%.2f |",sol);
-		} else if (p == 0)
-			sol = y - x;
-		else if (p == -0.5)
-			sol = asinh((y - x) / fabs(l));
-		else if (p == -1)
-			sol = atan((y - x) / fabs(l));
-		else if (p == -1.5)
-			sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l));
-		else
-			sol = (y - x) * pow((y - x) * (y - x) + l * l, p) + 2 * p * l * l
-					* g0(p - 1, y, x, l) / (2 * p + 1);
-	} else {
-		if (p == -0.5)
-			sol = sign(y - x) * log(fabs(y - x));
-		else
-			sol = (y - x) * pow(fabs(y - x), 2 * p) / (2 * p + 1);
-	}
-
-	return sol;
-}
-
-double G00(double p, double y1, double y2, double x1, double x2, double l) {
-	//	printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);
-
-	double sol = 0;
-	if (p == -1.5) {
-		if (l == 0) {
-			sol = -sqrt((y1 - x1) * (y1 - x1) + (y2 - x2) * (y2 - x2)) / ((y1
-					- x1) * (y2 - x2));
-		} else {
-			sol = sign((y1 - x1) * (y2 - x2)) / (2 * fabs(l)) * acos(
-					-2 * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2) / (((y1
-							- x1) * (y1 - x1) + l * l) * ((y2 - x2) * (y2 - x2)
-							+ l * l)) + 1);
-		}
-
-	} else if (p == -0.5) {
-		if (l != 0)
-			sol = 2 * p * l * l * G00(p - 1, y1, y2, x1, x2, l);
-		if ((y1 - x1) != 0)
-			sol += (y1 - x1) * g0(p, y2, x2,
-					sqrt((y1 - x1) * (y1 - x1) + l * l));
-		if ((y2 - x2) != 0)
-			sol += (y2 - x2) * g0(p, y1, x1,
-					sqrt((y2 - x2) * (y2 - x2) + l * l));
-		sol /= 2 * p + 2;
-	} else {
-		mexErrMsgTxt("no case for p defined");
-	}
-
-	return sol;
-}
-
-double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
-		double d3) {
-
-	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);
-	double sol = (x1 - y1 - d1) * (x2 - y2 - d2);
-
-	if (sol != 0)
-		sol *= slpADLO(x1, x2, y1 + d1, y2 + d2, d3);
-
-	if ((x1 - y1 - d1) != 0)
-		sol -= (x1 - y1 - d1) * compute_g0(0.5, x1, y1 + d1,
-				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));
-
-	if ((x2 - y2 - d2) != 0)
-		sol -= (x2 - y2 - d2) * compute_g0(0.5, x2, y2 + d2,
-				sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));
-
-	double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2
-			- d2) + d3 * d3);
-	sol += 1. / 3 * hlp * sqrt(hlp);
-	return sol;
-}
-
-double inline quadInt(
-		double(*f)(double, double, double, double, double, double, double),
-		double s1, double s2, double k1, double k2, double t1, double t2,
-		double l1, double l2, double d1, double d2, double d3) {
-
-	return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f(
-			k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1,
-			s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2,
-			t1, l2, d1, d2, d3) - f(k1, s2, t1, t2, d1, d2, d3) - f(s1, k2, l1,
-			l2, d1, d2, d3) + f(s1, k2, l1, t2, d1, d2, d3) + f(s1, k2, t1, l2,
-			d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1,
-			d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2,
-			d3) + f(s1, s2, t1, t2, d1, d2, d3);
-}
-
-double slpADLO(double y1, double y2, double x1, double x2, double a) {
-	double G3 = 0;
-	double gL = 0;
-	double gK = 0;
-	double tmp;
-
-	tmp = y1 - x1;
-	if (fabs(tmp) >= EPS * y1) {
-		tmp = sqrt(y1 * y1 + x1 * x1 + a * a - 2 * y1 * x1);
-		gL = compute_g0(-0.5, y2, x2, tmp);
-	}
-	tmp = y2 - x2;
-	if (fabs(tmp) >= EPS * y2) {
-		tmp = sqrt(y2 * y2 + x2 * x2 + a * a - 2 * y2 * x2);
-		gK = compute_g0(-0.5, y1, x1, tmp);
-	}
-	if (fabs(a * a) > EPS) {
-		if ((y1 - x1) * (y2 - x2) * a >= 0)
-			tmp = 1.;
-		else
-			tmp = -1.;
-
-		G3 = tmp * acos(
-				(-2. * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2)) / (((y1
-						- x1) * (y1 - x1) + a * a) * ((y2 - x2) * (y2 - x2) + a
-						* a)) + 1.) / (2. * a);
-	}
-
-	return (y1 - x1) * gL + (y2 - x2) * gK - a * a * G3;
-}
-
-double compute_g0(double p, double y, double x, double a) {
-	int sp = (int) 2 * (p - EPS); // MK (p-EPS) instead of (p+EPS)
-	//printf("\n compute_g0, p = %lf, sp = %d\n",p,sp);
-	assert(
-			p == 0 || (p == -0.5) || (p == -1) || (p == -1.5) || (fabs(a)
-					<= EPS));
-	if (fabs(a) <= EPS) {
-		// printf("\n a < eps\n");
-		switch (sp) {
-		case 0:
-			return y - x;
-		case -1:
-			return log(fabs(x - y)) * (y - x) / fabs(y - x);
-		case -2:
-			return -(y - x) / fabs(y - x) / fabs(y - x);
-		case -3:
-			return -0.5 * (y - x) / fabs(y - x) / fabs(y - x) / fabs(y - x);
-		}
-	} else {
-		//  printf("\n a > eps\n");
-		switch (sp) {
-		case 0:
-			return y - x;
-		case -1:
-			return asinh((y - x) / fabs(a));
-		case -2:
-			return atan((y - x) / fabs(a));
-		case -3:
-			return (y - x) * pow((x * x + y * y + a * a - 2 * x * y), -0.5)
-					/ (a * a);
-		default:
-			printf("p must be either 0, -1/2, -1 or -3/2.");
-			return NAN;
-		}
-	}
-}
-
-double FLO_plane(double x1, double x2, double y1, double y2, double delta1,
-		double delta2, double a) {
-	double yd1 = y1 + delta1;
-	double yd2 = y2 + delta2;
-	double tmp1 = x1 - y1 - delta1;
-	double tmp2 = x2 - y2 - delta2;
-	double tmp3 = sqrt(tmp2 * tmp2 + a * a);
-	double tmp4 = sqrt(tmp1 * tmp1 + a * a);
-	double tmp5 = pow(tmp1 * tmp1 + tmp2 * tmp2 + a * a, 3. / 2.);
-	double rval = 0;
-
-	rval = tmp1 * tmp2 * slpADLO(x1, x2, yd1, yd2, a) - tmp1 * compute_g0(-0.5,
-			x1, yd1, tmp3) - tmp2 * compute_g0(-0.5, x2, yd2, tmp4) + tmp5 / 3.;
-	return rval;
-}
diff --git a/src/bem3d/src/laplace/SLPrecangle.hpp b/src/bem3d/src/laplace/SLPrecangle.hpp
deleted file mode 100644
index 6f80150..0000000
--- a/src/bem3d/src/laplace/SLPrecangle.hpp
+++ /dev/null
@@ -1,21 +0,0 @@
-#ifndef HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
-#define HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
-
-int sign(double);
-double arsinh(double);
-
-// sol = g0(p,y,x,l);
-double g0(double, double, double, double);
-// sol = G00(p,y1,y2,x1,x2,l);
-double G00(double, double, double, double, double, double);
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
-double F_par(double, double, double, double, double, double, double);
-// sol = F_ort(x1,x2,y1,y2,d1,d2,d3);
-//double F_ort(double, double, double, double, double, double, double);
-
-
-double slpADLO(double, double, double, double, double);
-double compute_g0(double, double, double, double);
-double FLO_plane(double, double, double, double, double, double, double);
-
-#endif
diff --git a/src/bem3d/src/laplace/SLPrectangle.cpp b/src/bem3d/src/laplace/SLPrectangle.cpp
new file mode 100644
index 0000000..26db04b
--- /dev/null
+++ b/src/bem3d/src/laplace/SLPrectangle.cpp
@@ -0,0 +1,166 @@
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include <stdlib.h>
+
+#include "SLPrectangle.hpp"
+
+#define my_PI 3.141592653589793
+#define EPS 0.00001
+
+using namespace std;
+
+int sign(double);
+//double arsinh(double);
+
+
+int inline sign(double x) {
+	return x > 0 ? 1 : (x < 0 ? -1 : 0);
+}
+/*
+ double inline arsinh(double x) {
+ return log(x + sqrt(x * x + 1));
+ }
+ */
+
+//y-x muss != 0 sein
+double g0(double p, double y, double x, double l) {
+	//printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);
+
+	double sol = 0;
+
+	if (l != 0) {
+		if (p == 0.5) {
+			sol = (y - x) / 2 * sqrt((y - x) * (y - x) + l * l) + l * l / 2
+					* asinh((y - x) / fabs(l));
+			//          printf("%.2f |",sol);
+		} else if (p == 0)
+			sol = y - x;
+		else if (p == -0.5)
+			sol = asinh((y - x) / fabs(l));
+		else if (p == -1)
+			sol = atan((y - x) / fabs(l)) / l;
+		else if (p == -1.5)
+			sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l));
+		else
+			sol = (y - x) * pow((y - x) * (y - x) + l * l, p) + 2 * p * l * l
+					* g0(p - 1, y, x, l) / (2 * p + 1);
+	} else {
+		if (p == -0.5)
+			sol = sign(y - x) * log(fabs(y - x));
+		else
+			sol = (y - x) * pow(fabs(y - x), 2 * p) / (2 * p + 1);
+	}
+
+	return sol;
+}
+
+double G00(double p, double y1, double y2, double x1, double x2, double l) {
+	//	printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);
+
+	double sol = 0;
+	if (p == -1.5) {
+		if (l == 0) {
+			sol = -sqrt((y1 - x1) * (y1 - x1) + (y2 - x2) * (y2 - x2)) / ((y1
+					- x1) * (y2 - x2));
+		} else {
+			sol = sign((y1 - x1) * (y2 - x2)) / (2 * fabs(l)) * acos(
+					-2 * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2) / (((y1
+							- x1) * (y1 - x1) + l * l) * ((y2 - x2) * (y2 - x2)
+							+ l * l)) + 1);
+		}
+
+	} else if ((pt > 0) && ((int) pt == pt)) {
+		if (l != 0)
+			sol = 2 * p * l * l * G00(p - 1, y1, y2, x1, x2, l);
+		if ((y1 - x1) != 0)
+			sol += (y1 - x1) * g0(p, y2, x2,
+					sqrt((y1 - x1) * (y1 - x1) + l * l));
+		if ((y2 - x2) != 0)
+			sol += (y2 - x2) * g0(p, y1, x1,
+					sqrt((y2 - x2) * (y2 - x2) + l * l));
+		sol /= 2 * p + 2;
+	} else {
+		sol = NAN;
+		printf("warning in G00: no case for p=%.2f defined\n", p);
+	}
+
+	return sol;
+}
+
+double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
+		double d3) {
+
+	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);
+	double sol = (x1 - y1 - d1) * (x2 - y2 - d2);
+
+	if (sol != 0)
+		sol *= G00(-0.5, x1, x2, y1 + d1, y2 + d2, d3);
+
+	if ((x1 - y1 - d1) != 0)
+		sol -= (x1 - y1 - d1) * g0(0.5, x1, y1 + d1,
+				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));
+
+	if ((x2 - y2 - d2) != 0)
+		sol -= (x2 - y2 - d2) * g0(0.5, x2, y2 + d2,
+				sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));
+
+	double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2
+			- d2) + d3 * d3);
+	sol += 1. / 3 * hlp * sqrt(hlp);
+	return sol;
+}
+
+double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,
+		double d3) {
+
+	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y2,y3,d1,d2,d3);
+	double sol = -G00(0.5, y3, x1, -d3, d1, x2 - y2 - d2);
+
+	if ((x1 - d1) * (x2 - y2 - d2) != 0)
+		sol -= (x1 - d1) * (x2 - y2 - d2) * G00(-0.5, x2, y3, y2 + d2, -d3,
+				x1 - d1);
+
+	if ((x1 - d1) != 0)
+		sol += (x1 - d1) * g0(0.5, y3, -d3,
+				sqrt((x1 - d1) * (x1 - d1) + (x2 - y2 - d2) * (x2 - y2 - d2)));
+
+	if ((y3 + d3) * (x2 - y2 - d2) != 0)
+		sol -= (y3 + d3) * (x2 - y2 - d2) * G00(-0.5, x1, x2, d1, y2 + d2,
+				-y3 - d3);
+
+	if ((y3 + d3) != 0)
+		sol += (y3 + d3) * g0(0.5, x1, d1,
+				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + (y3 + d3) * (y3 + d3)));
+
+	return sol / 2.;
+}
+
+double quadInt(
+		double(*f)(double, double, double, double, double, double, double),
+		double a, double b, double c, double d, double r, double t, double u,
+		double v, double d1, double d2, double d3) {
+
+	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) - 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) - 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)
+			+ 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);
+}
+
+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/bem3d/src/laplace/SLPrectangle.hpp b/src/bem3d/src/laplace/SLPrectangle.hpp
new file mode 100644
index 0000000..eaff0f0
--- /dev/null
+++ b/src/bem3d/src/laplace/SLPrectangle.hpp
@@ -0,0 +1,28 @@
+#ifndef HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
+#define HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
+
+//int sign(double);
+//double arsinh(double);
+
+// sol = g0(p,y,x,l);
+double g0(double, double, double, double);
+// sol = G00(p,y1,y2,x1,x2,l);
+double G00(double, double, double, double, double, double);
+// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
+double F_par(double, double, double, double, double, double, double);
+// sol = F_ort(x1,x2,y2,y3,d1,d2,d3);
+double F_ort(double, double, double, double, double, double, double);
+
+// sol = quadInt((F_par/F_ort),a,b,c,d,r,t,u,v,d1,d2,d3);
+double quadInt(
+		double(*)(double, double, double, double, double, double, double),
+		double, double, double, double, double, double, double, double, double,
+		double, double);
+
+// sol = quad0Int((F_par/F_ort),b,d,t,v,d1,d2,d3);
+double quad0Int(
+		double(*)(double, double, double, double, double, double, double),
+		double, double, double, double, double, double, double);
+
+
+#endif
diff --git a/src/bem3d/src/laplace/SimpleLayerPotential.cpp b/src/bem3d/src/laplace/SimpleLayerPotential.cpp
index c5e299a..8fd3573 100644
--- a/src/bem3d/src/laplace/SimpleLayerPotential.cpp
+++ b/src/bem3d/src/laplace/SimpleLayerPotential.cpp
@@ -3,174 +3,284 @@
 #include "../spaces/P0Space.hpp"
 #include "TriangleIntegrator.hpp"
 
-#include "SLPrecangle.hpp"
-
+#include "SLPrectangle.hpp"	// Bibliothek fuer achsenorientierte Rechtecke
 #include <fstream>
 
 using namespace laplace;
 
-SimpleLayerPotential::SimpleLayerPotential(const space::P0Space& space)
-  : HLibOperator< space::P0Space, space::P0Space, SymmetricMatrix >(
-        space, space, SimpleLayerPotential::computeEntry)
-{
+SimpleLayerPotential::SimpleLayerPotential(const space::P0Space& space) :
+			HLibOperator<space::P0Space, space::P0Space, SymmetricMatrix> (
+					space, space, SimpleLayerPotential::computeEntry) {
 }
 
-SimpleLayerPotential::~SimpleLayerPotential()
-{
+SimpleLayerPotential::~SimpleLayerPotential() {
 }
 
-template < class TMesh >
-double ComputeSimpleLayerPotentialEntry< TMesh >::exec(
-    const typename TMesh::Prop::Triangle& f1,
-    const typename TMesh::Prop::Triangle& f2) const
-{
-  TriangleIntegrator integrator = TriangleIntegrator();
-
-  if ( f1.computeArea() < f2.computeArea() )
-  {
-    return integrator.integrate( f1, f2,
-          computeInnerIntegralForTriangles);
-  }
-  else
-  {
-    return integrator.integrate( f2, f1,
-          computeInnerIntegralForTriangles);
-  }
+template<class TMesh>
+double ComputeSimpleLayerPotentialEntry<TMesh>::exec(
+		const typename TMesh::Prop::Triangle& f1,
+		const typename TMesh::Prop::Triangle& f2) const {
+	TriangleIntegrator integrator = TriangleIntegrator();
+
+	if (f1.computeArea() < f2.computeArea()) {
+		return integrator.integrate(f1, f2, computeInnerIntegralForTriangles);
+	} else {
+		return integrator.integrate(f2, f1, computeInnerIntegralForTriangles);
+	}
 }
 
-template < class TMesh >
-double ComputeSimpleLayerPotentialEntry< TMesh >::exec(
-    const typename TMesh::Prop::Parallelogram& f1,
-    const typename TMesh::Prop::Parallelogram& f2) const
-{
-
-   const typename TMesh::Vertex& a1 = f1.getA();
-   const typename TMesh::Vertex& a2 = f2.getA();
-
-
-
-
-   return NAN;
+template<class TMesh>
+double ComputeSimpleLayerPotentialEntry<TMesh>::exec(
+		const typename TMesh::Prop::Parallelogram& f1,
+		const typename TMesh::Prop::Parallelogram& f2) const {
+
+	// Zwischenvariable fuer die Loesung (kann fuer geschwindigkeit entfernt werden)
+	double sol = NAN;
+
+	// Flaechenstueck 1
+	const typename TMesh::Vertex& xa = f1.getA() - f1.getB();
+	const typename TMesh::Vertex& xb = f1.getA() - f1.getD();
+	const typename TMesh::Vertex& xn = f1.getA() * f1.getD(); //todo Kreuzprodukt!?
+
+	// Flaechenstueck 2
+	const typename TMesh::Vertex& ya = f2.getA() - f1.getB();
+	const typename TMesh::Vertex& yb = f2.getA() - f1.getD();
+	const typename TMesh::Vertex& yn = f2.getA() * f2.getD(); //todo Kreuzprodukt!?
+
+	// Abstand zwischen den Stuecken (Achtung) abstand von den beiden "kleinsten" Ecken
+	const typename TMesh::Vertex& d = NULL; //todo new???
+
+	//Variablen fuer LageInformationen
+	int rx, rxa, rxb, ry, rya, ryb;
+
+	// Lage des ersten Elements ermitteln
+	if (xn[2] != 0)
+		rx = 2;
+	else if (xn[1] != 0)
+		rx = 1;
+	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;
+
+	//kleinste Ecke finden und fuer \delta verwenden
+	if (xa[rxa] > 0) {
+		if (xb[rxb] > 0) {
+			d = -f1.getA;
+		} else {
+			d = -f1.getD;
+		}
+	} else {
+		if (xb[rxb] > 0) {
+			d = -f1.getB;
+		} else {
+			d = -f1.getC;
+		}
+	}
+
+	// Lage des zweiten Elements ermitteln
+	if (yn[2] != 0)
+		ry = 2;
+	else if (yn[1] != 0)
+		ry = 1;
+	else
+		ry = 0;
+
+	if (ya[2] != 0)
+		rya = 2;
+	else if (ya[1] != 0)
+		rya = 1;
+	else
+		rya = 0;
+
+	if (yb[2] != 0)
+		ryb = 2;
+	else if (yb[1] != 0)
+		ryb = 1;
+	else
+		ryb = 0;
+
+	//kleinste Ecke finden und fuer \delta verwenden
+	if (ya[rya] > 0) {
+		if (yb[ryb] > 0) {
+			d += f2.getA; //todo += definiert??? sonst d = d + ...
+		} else {
+			d += f2.getD; //todo += definiert??? sonst d = d + ...
+		}
+	} else {
+		if (yb[ryb] > 0) {
+			d += f2.getB; //todo += definiert??? sonst d = d + ...
+		} else {
+			d += f2.getC; //todo += definiert??? sonst d = d + ...
+		}
+	}
+
+	if (rx == ry) { // Flaechen sind parallel
+		if (rxa == rya) { // Elemente zeigen in gleiche Richtung
+			sol = quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]),
+					fabs(yb[rxb]), d[rxa], d[rxb], d[rx]);
+		} else { // Elemente zeigen nicht in gleiche Richtung
+			sol = quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[rxa]),
+					fabs(ya[rxb]), d[rxa], d[rxb], d[rx]);
+		}
+	} else { // Flaechen sind orthogonal
+		if (rxa == rya) { //Elmente richtig zuordnen
+			sol = quad0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]), fabs(ya[rya]),
+					fabs(yb[ryb]), d[rxb], d[rxa], d[rx]);
+		} else if (rxa == ryb) {
+			sol = quad0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]), fabs(yb[ryb]),
+					fabs(ya[rya]), d[rxb], d[rxa], d[rx]);
+		} else if (rxb == rya) {
+			sol = quad0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rya]),
+					fabs(yb[ryb]), d[rxa], d[rxb], d[rx]);
+		} else {
+			sol = quad0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[ryb]),
+					fabs(ya[rya]), d[rxa], d[rxb], d[rx]);
+		}
+		sol /= 2.;
+	}
+
+	return sol / (4. * PI);
 }
 
-double SimpleLayerPotential::computeEntry(
-    const space::P0Space::Element &e1,
-    const space::P0Space::Element &e2)
-{
-  typedef boundary_mesh::DoubleFaceDispatcher<
-    ComputeSimpleLayerPotentialEntry< BoundaryMesh >, BoundaryMesh, double >
-  ComputeEntry;
+double SimpleLayerPotential::computeEntry(const space::P0Space::Element &e1,
+		const space::P0Space::Element &e2) {
+	typedef boundary_mesh::DoubleFaceDispatcher<
+			ComputeSimpleLayerPotentialEntry<BoundaryMesh> , BoundaryMesh,
+			double> ComputeEntry;
 
-  return ComputeEntry::Exec( e1.getSupport(), e2.getSupport(),
-                      ComputeSimpleLayerPotentialEntry< BoundaryMesh >() );
+	return ComputeEntry::Exec(e1.getSupport(), e2.getSupport(),
+			ComputeSimpleLayerPotentialEntry<BoundaryMesh> ());
 }
 
-template < class TMesh >
-double
-ComputeSimpleLayerPotentialEntry< TMesh >::computeInnerIntegralForTriangles(
-    const typename TMesh::Prop::Triangle& tau,
-    const typename TMesh::UVertex& x)
-{
-  BoundaryMesh::Face::ConstVertexIterator vIt = tau.begin_v();
-  const BoundaryMesh::UVertex& x1 = *vIt; ++vIt;
-  const BoundaryMesh::UVertex& x2 = *vIt; ++vIt;
-  const BoundaryMesh::UVertex& x3 = *vIt; ++vIt;
-
-  // determine a local orthogonal coordinate system (r1,r2,n) on the triangle tau
-  double t_tau = (x3-x2).norm();
-  BoundaryMesh::UVertex r2 = (x3-x2)/t_tau;
-  double t_star=(x1-x2)*r2;
-
-  BoundaryMesh::UVertex r1 = (x2+t_star*r2-x1);
-  double s_tau = r1.norm();
-  r1 = r1 / s_tau;
-
-  BoundaryMesh::UVertex n = BoundaryMesh::Vertex::crossProduct(r1,r2);
-
-  // determine necessary information
-  double sx = (x-x1)*r1;
-  double tx = (x-x1)*r2;
-  double ux = (x-x1)*n;
-  double a1 = -t_star/s_tau;
-  double a2 = (t_tau-t_star)/s_tau;
-
-  // analytic computation of inner integral
-  return 0.25*(F(s_tau,a2,sx,tx,ux) - F(0,a2,sx,tx,ux) - F(s_tau,a1,sx,tx,ux) + F(0,a1,sx,tx,ux))/M_PI;
+template<class TMesh>
+double ComputeSimpleLayerPotentialEntry<TMesh>::computeInnerIntegralForTriangles(
+		const typename TMesh::Prop::Triangle& tau,
+		const typename TMesh::UVertex& x) {
+	BoundaryMesh::Face::ConstVertexIterator vIt = tau.begin_v();
+	const BoundaryMesh::UVertex& x1 = *vIt;
+	++vIt;
+	const BoundaryMesh::UVertex& x2 = *vIt;
+	++vIt;
+	const BoundaryMesh::UVertex& x3 = *vIt;
+	++vIt;
+
+	// determine a local orthogonal coordinate system (r1,r2,n) on the triangle tau
+	double t_tau = (x3 - x2).norm();
+	BoundaryMesh::UVertex r2 = (x3 - x2) / t_tau;
+	double t_star = (x1 - x2) * r2;
+
+	BoundaryMesh::UVertex r1 = (x2 + t_star * r2 - x1);
+	double s_tau = r1.norm();
+	r1 = r1 / s_tau;
+
+	BoundaryMesh::UVertex n = BoundaryMesh::Vertex::crossProduct(r1, r2);
+
+	// determine necessary information
+	double sx = (x - x1) * r1;
+	double tx = (x - x1) * r2;
+	double ux = (x - x1) * n;
+	double a1 = -t_star / s_tau;
+	double a2 = (t_tau - t_star) / s_tau;
+
+	// analytic computation of inner integral
+	return 0.25 * (F(s_tau, a2, sx, tx, ux) - F(0, a2, sx, tx, ux) - F(s_tau,
+			a1, sx, tx, ux) + F(0, a1, sx, tx, ux)) / M_PI;
 }
 
-template < class TMesh >
-double
-ComputeSimpleLayerPotentialEntry< TMesh >::F(
-  double s, double a, double sx, double tx, double ux)
-{
-  double p = (a*tx+sx)/(1+a*a);
-  double q = ux*ux + (tx-a*sx)*(tx-a*sx)/(1+a*a);
-
-  return (s-sx)*log(a*s - tx + sqrt((s-sx)*(s-sx) + (a*s-tx)*(a*s-tx) + ux*ux)) - s
-    + (a*sx-tx)*log(sqrt(1+a*a)*(s-p)+sqrt((1+a*a)*(s-p)*(s-p) + q*q))/sqrt(q+a*a)
-    + 2*ux*atan(((q-(a*sx-tx)/(1+a*a))*sqrt((1+a*a)*(s-p)*(s-p)+q*q) + (a*s-tx-q)*q)/((s-p)*ux));
+template<class TMesh>
+double ComputeSimpleLayerPotentialEntry<TMesh>::F(double s, double a,
+		double sx, double tx, double ux) {
+	double p = (a * tx + sx) / (1 + a * a);
+	double q = ux * ux + (tx - a * sx) * (tx - a * sx) / (1 + a * a);
+
+	return (s - sx) * log(
+			a * s - tx
+					+ sqrt(
+							(s - sx) * (s - sx) + (a * s - tx) * (a * s - tx)
+									+ ux * ux)) - s + (a * sx - tx) * log(
+			sqrt(1 + a * a) * (s - p) + sqrt(
+					(1 + a * a) * (s - p) * (s - p) + q * q)) / sqrt(q + a * a)
+			+ 2 * ux * atan(
+					((q - (a * sx - tx) / (1 + a * a)) * sqrt(
+							(1 + a * a) * (s - p) * (s - p) + q * q) + (a * s
+							- tx - q) * q) / ((s - p) * ux));
 }
 
-void fill_supermatrix_aca_OLD(supermatrix* A, const cluster* row, const cluster* col, void *data,
-    double (*compute_entry)(int i, int j, void* data), double aca_eps)
-{
-//  std::cerr << "*** fill_supermatrix_aca" << std::endl;
-
-  fullmatrix* f = A->f;
-  rkmatrix* r = A->r;
-  int i, j;
-  int rank;
-
-  if ( f!= NULL ) { // fill fullmatrix block
-//    std::cerr << "*** full matrix block" << std::endl;
-    for ( i=0; i<row->size; i++ ){
-      for ( j=0; j<col->size; j++ )
-	f->e[i+j*row->size] = compute_entry(row->start+i,col->start+j, data);
-    }
-
-    return;
-  }
-
-  if ( r!= NULL ) { // fill rk-block with aca
-//    std::cerr << "*** rank k matrix block" << std::endl;
-    if ( r->rows > r-> cols )
-      rank = r->cols;
-    else
-      rank = r->rows;
-
-    if ( rank > ACA_KMAX )
-      rank = ACA_KMAX;
-    
-    r->k = rank;
-    r->a = allocate_matrix(r->rows,rank);
-    r->b = allocate_matrix(r->cols,rank);
-	
-    r->kt = newaca_fill_block(r->a, r->b, r->rows, r->cols,
-				  row->start, col->start, 
-				  compute_entry, data,
-				  rank, aca_eps, ACA_STRATEGY);
-
-    r->a = (double*) realloc(r->a,r->rows*r->kt*sizeof(double));
-    r->b = (double*) realloc(r->b,r->cols*r->kt*sizeof(double));
-    r->k = r->kt;
-
-    return;
-  }
-
-  if(row->sons>0) {
-    if(col->sons>0) {
-      for(int i=0;i<row->sons;++i)
-        for(int j=0;j<col->sons;++j) {
-            int s=i+j*row->sons;
-            fill_supermatrix_aca_OLD(A->s[s],row->son[i],col->son[j],data,compute_entry,aca_eps);
-        }
-    } else {
-      for(int i=0;i<row->sons;++i)
-        fill_supermatrix_aca_OLD(A->s[i],row->son[i],col,data,compute_entry,aca_eps);
-    }
-  } else {
-    for(int i=0;i<col->sons;++i)
-      fill_supermatrix_aca_OLD(A->s[i],row,col->son[i],data,compute_entry,aca_eps);
-  }
+void fill_supermatrix_aca_OLD(supermatrix* A, const cluster* row,
+		const cluster* col, void *data,
+		double(*compute_entry)(int i, int j, void* data), double aca_eps) {
+	//  std::cerr << "*** fill_supermatrix_aca" << std::endl;
+
+	fullmatrix* f = A->f;
+	rkmatrix* r = A->r;
+	int i, j;
+	int rank;
+
+	if (f != NULL) { // fill fullmatrix block
+	//    std::cerr << "*** full matrix block" << std::endl;
+		for (i = 0; i < row->size; i++) {
+			for (j = 0; j < col->size; j++)
+				f->e[i + j * row->size] = compute_entry(row->start + i,
+						col->start + j, data);
+		}
+
+		return;
+	}
+
+	if (r != NULL) { // fill rk-block with aca
+	//    std::cerr << "*** rank k matrix block" << std::endl;
+		if (r->rows > r-> cols)
+			rank = r->cols;
+		else
+			rank = r->rows;
+
+		if (rank > ACA_KMAX)
+			rank = ACA_KMAX;
+
+		r->k = rank;
+		r->a = allocate_matrix(r->rows, rank);
+		r->b = allocate_matrix(r->cols, rank);
+
+		r->kt = newaca_fill_block(r->a, r->b, r->rows, r->cols, row->start,
+				col->start, compute_entry, data, rank, aca_eps, ACA_STRATEGY);
+
+		r->a = (double*) realloc(r->a, r->rows * r->kt * sizeof(double));
+		r->b = (double*) realloc(r->b, r->cols * r->kt * sizeof(double));
+		r->k = r->kt;
+
+		return;
+	}
+
+	if (row->sons > 0) {
+		if (col->sons > 0) {
+			for (int i = 0; i < row->sons; ++i)
+				for (int j = 0; j < col->sons; ++j) {
+					int s = i + j * row->sons;
+					fill_supermatrix_aca_OLD(A->s[s], row->son[i], col->son[j],
+							data, compute_entry, aca_eps);
+				}
+		} else {
+			for (int i = 0; i < row->sons; ++i)
+				fill_supermatrix_aca_OLD(A->s[i], row->son[i], col, data,
+						compute_entry, aca_eps);
+		}
+	} else {
+		for (int i = 0; i < col->sons; ++i)
+			fill_supermatrix_aca_OLD(A->s[i], row, col->son[i], data,
+					compute_entry, aca_eps);
+	}
 }
 
diff --git a/src/bem3d/src/spaces/testP0Space b/src/bem3d/src/spaces/testP0Space
index fc3fe2f..09b9410 100755
Binary files a/src/bem3d/src/spaces/testP0Space and b/src/bem3d/src/spaces/testP0Space differ
diff --git a/src/build_A.cpp b/src/build_A.cpp
index 4bf70d2..91d5c82 100644
--- a/src/build_A.cpp
+++ b/src/build_A.cpp
@@ -16,7 +16,7 @@ using namespace std;
 
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
-	int i,j,k;
+	int i, j, k;
 
 	//sicherheitsabfragen zu Datengroessen
 	if (nrhs != 2)
@@ -82,10 +82,10 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 		x3[1] = C[cm + (int) E[3 * em + j] - 1];
 		x3[2] = C[2 * cm + (int) E[3 * em + j] - 1];
 
-		for (i = 0; i<3;++i)
+		for (i = 0; i < 3; ++i)
 			xa[i] = x1[i] - x0[i];
 
-		for (i = 0; i<3;++i)
+		for (i = 0; i < 3; ++i)
 			xb[i] = x3[i] - x0[i];
 
 		//Kreuzprodukt
@@ -93,7 +93,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 		xn[1] = (xa[2]) * (xb[0]) - (xa[0]) * (xb[2]);
 		xn[2] = (xa[0]) * (xb[1]) - (xa[1]) * (xb[0]);
 
-
 		// Lageeigenschaften des Flaechenstuecks
 		//if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])
 		if (xn[2] != 0)
@@ -119,21 +118,21 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 		//kleinste Ecke finden und fuer \delta verwenden
 
-		if(xa[rxa]>0){
-			if(xb[rxb]>0){
-				for (i = 0; i<3;++i)
-					dt[i] = - x0[i];
-			}else{
-				for (i = 0; i<3;++i)
-					dt[i] = - x3[i];
+		if (xa[rxa] > 0) {
+			if (xb[rxb] > 0) {
+				for (i = 0; i < 3; ++i)
+					dt[i] = -x0[i];
+			} else {
+				for (i = 0; i < 3; ++i)
+					dt[i] = -x3[i];
 			}
-		}else{
-			if(xb[rxb]>0){
-				for (i = 0; i<3;++i)
-					dt[i] = - x1[i];
-			}else{
-				for (i = 0; i<3;++i)
-					dt[i] = - x2[i];
+		} else {
+			if (xb[rxb] > 0) {
+				for (i = 0; i < 3; ++i)
+					dt[i] = -x1[i];
+			} else {
+				for (i = 0; i < 3; ++i)
+					dt[i] = -x2[i];
 				//printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");
 			}
 		}
@@ -161,10 +160,10 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 			y3[1] = C[cm + (int) E[3 * em + k] - 1];
 			y3[2] = C[2 * cm + (int) E[3 * em + k] - 1];
 
-			for (i = 0; i<3;++i)
+			for (i = 0; i < 3; ++i)
 				ya[i] = y1[i] - y0[i];
 
-			for (i = 0; i<3;++i)
+			for (i = 0; i < 3; ++i)
 				yb[i] = y3[i] - y0[i];
 
 			yn[0] = (ya[1]) * (yb[2]) - (ya[2]) * (yb[1]);
@@ -194,27 +193,27 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 			//kleinste Ecke finden und fuer \delta verwenden
 
-			if(ya[rya]>0){
-				if(yb[ryb]>0){
-					for (i = 0; i<3;++i)
+			if (ya[rya] > 0) {
+				if (yb[ryb] > 0) {
+					for (i = 0; i < 3; ++i)
 						d[i] = dt[i] + y0[i];
-				}else{
-					for (i = 0; i<3;++i)
+				} else {
+					for (i = 0; i < 3; ++i)
 						d[i] = dt[i] + y3[i];
 				}
-			}else{
-				if(yb[ryb]>0){
-					for (i = 0; i<3;++i)
+			} else {
+				if (yb[ryb] > 0) {
+					for (i = 0; i < 3; ++i)
 						d[i] = dt[i] + y1[i];
-				}else{
-					for (i = 0; i<3;++i)
+				} else {
+					for (i = 0; i < 3; ++i)
 						d[i] = dt[i] + y2[i];
 					//printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");
 				}
 			}
 
-		//	for (i = 0; i<3;++i)
-		//		d[i] = y0[i]-x0[i];
+			//	for (i = 0; i<3;++i)
+			//		d[i] = y0[i]-x0[i];
 
 			//printf("(%d,%d)",rx,ry);
 			if (rx == ry) {
@@ -226,22 +225,52 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],
 					//		y1[0], y0[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 = quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),
-							fabs(ya[rxa]), fabs(yb[rxb]), d[rxa], d[rxb], d[rx]);
+					tmp
+							= quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(ya[rxa]), fabs(yb[rxb]), d[rxa],
+									d[rxb], d[rx]);
 					//	 printf("%d%d|%.2f\n",j,k,tmp);
-				} else{
+				} else {
 					//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],
 					//		y1[0], y0[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 = quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),
-							fabs(yb[rxa]), fabs(ya[rxb]), d[rxa], d[rxb], d[rx]);
+					tmp
+							= quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(yb[rxa]), fabs(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;
 
 			} else {
-				A[(k * em) + j] = NAN;
+				//printf("(%d,%d@%d,%d)", rxa, rxb, rya, ryb);
+				//printf("(%d,%d)", rx, ry);
+				//printf("\n");
+
+				if (rxa == rya) {
+					tmp
+							= quad0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),
+									fabs(ya[rya]), fabs(yb[ryb]), d[rxb],
+									d[rxa], d[rx]);
+				} else if (rxa == ryb) {
+					tmp
+							= quad0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),
+									fabs(yb[ryb]), fabs(ya[rya]), d[rxb],
+									d[rxa], d[rx]);
+				} else if (rxb == rya) {
+					tmp
+							= quad0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(ya[rya]), fabs(yb[ryb]), d[rxa],
+									d[rxb], d[rx]);
+				} else {
+					tmp
+							= quad0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(yb[ryb]), fabs(ya[rya]), d[rxa],
+									d[rxb], d[rx]);
+				}
+				A[(k * em) + j] = tmp;
+
 			}
+			A[(k * em) + j] = 1. / (4 * my_PI) * tmp;
 
 			// Vorbereiten der DATEN
 			//   ej = coordinates(elements(j,:)',:);
diff --git a/src/exmpl_3DFichCube.mat b/src/exmpl_3DFichCube.mat
new file mode 100644
index 0000000..3c5c502
Binary files /dev/null and b/src/exmpl_3DFichCube.mat differ
diff --git a/src/exportCOEL.m b/src/exportCOEL.m
index 7dbe086..30fe941 100644
--- a/src/exportCOEL.m
+++ b/src/exportCOEL.m
@@ -29,18 +29,18 @@ if(size(ele,2) == 3)
   end
 end
 if(size(ele,2)==4)
-%   disp 'schreibe in Datei'
-%   dat = fopen(file,'w');
-% 
-%       for i=1:size(coo,1);
-%           fprintf(dat,'v %.4f %.4f %.4f\n',coo(i,:));
-%       end
-%           fprintf(dat,'\n');
-%       for i=1:size(ele,1);
-%           fprintf(dat,'f %.0f %.0f %.0f %.0f\n',ele(i,:));
-%       end
-% 
-%   fclose(dat);
+  disp 'schreibe in Datei'
+  dat = fopen(file,'w');
+
+      for i=1:size(coo,1);
+          fprintf(dat,'v %.4f %.4f %.4f\n',coo(i,:));
+      end
+          fprintf(dat,'\n');
+      for i=1:size(ele,1);
+          fprintf(dat,'f %.0f %.0f %.0f %.0f\n',ele(i,:));
+      end
+
+  fclose(dat);
 end
 end
 
diff --git a/src/mex_Fort.cpp b/src/mex_Fort.cpp
new file mode 100644
index 0000000..7d4f894
--- /dev/null
+++ b/src/mex_Fort.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 (x1,x2,y2,y3,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_ort(*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_perp(*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_Fpar.cpp b/src/mex_Fpar.cpp
index 9af5fac..ab8ae94 100644
--- a/src/mex_Fpar.cpp
+++ b/src/mex_Fpar.cpp
@@ -15,7 +15,7 @@ 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)");
+		mexErrMsgTxt("expected (x1,x2,y1,y2,d1,d2,d3)");
 	if (nlhs > 2)
 		mexErrMsgTxt("has maximal two output arguments");
 
diff --git a/src/mexit.m b/src/mexit.m
index 4057572..f552eff 100644
--- a/src/mexit.m
+++ b/src/mexit.m
@@ -2,6 +2,7 @@
 mex mex_g0.cpp SLPrecangle.cpp
 mex mex_G00.cpp SLPrecangle.cpp
 mex mex_Fpar.cpp SLPrecangle.cpp
+mex mex_Fort.cpp SLPrecangle.cpp
 
 mex build_A.cpp SLPrecangle.cpp
 
diff --git a/src/surfDoubleQuad.m b/src/surfDoubleQuad.m
index a93accd..86ed43c 100644
--- a/src/surfDoubleQuad.m
+++ b/src/surfDoubleQuad.m
@@ -45,60 +45,60 @@
 % end
 
 
-% function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
-% 
-%     [x1n x1w] = gauss(w,a,b);
-%     [x2n x2w] = gauss(w,c,d);
-%     [y1n y1w] = gauss(w,r,t);
-%     [y2n y2w] = gauss(w,u,v);
-% 
-%     in = 0;
-%     for i=1:length(y2n)
-%         for j=1:length(y1n)
-%             for k=1:length(x2n)
-%                 for l = 1:length(x1n)
-%                     
-%                     in = in + y2w(i) * y1w(j) * x2w(k) * x1w(l) * f(x1n(l),x2n(k),y1n(j),y2n(i));
-%                     
-%                 end
-%             end
-%         end
-%     end
-% 
-% end
-
 function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
 
     [x1n x1w] = gauss(w,a,b);
     [x2n x2w] = gauss(w,c,d);
-    
-   
-     d1=r-a;
-     d2=u-c;
-%      sum(x1w)a
-     x1n=x1n-d1; x2n=x2n-d2;
-%     c
-%     d
-    
-    p = -0.5;
+    [y1n y1w] = gauss(w,r,t);
+    [y2n y2w] = gauss(w,u,v);
 
     in = 0;
-
+    for i=1:length(y2n)
+        for j=1:length(y1n)
             for k=1:length(x2n)
                 for l = 1:length(x1n)
                     
-    [my(1) mk(1)] = mex_G00(p,r,u,x1n(l),x2n(k),0);
-    [my(2) mk(2)] = mex_G00(p,t,v,x1n(l),x2n(k),0);
-    [my(3) mk(3)] = mex_G00(p,t,u,x1n(l),x2n(k),0);
-    [my(4) mk(4)] = mex_G00(p,r,v,x1n(l),x2n(k),0);
-
-
-%     my_sol = my(1)+my(2)-my(3)-my(4);
-                    
-                    in = in + x2w(k) * x1w(l) * (my(1)+my(2)-my(3)-my(4));
+                    in = in + y2w(i) * y1w(j) * x2w(k) * x1w(l) * f(x1n(l),x2n(k),y1n(j),y2n(i));
                     
                 end
             end
-%             in
+        end
+    end
+
+end
 
-end
\ No newline at end of file
+% function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
+% 
+%     [x1n x1w] = gauss(w,a,b);
+%     [x2n x2w] = gauss(w,c,d);
+%     
+%    
+%      d1=r-a;
+%      d2=u-c;
+% %      sum(x1w)a
+%      x1n=x1n-d1; x2n=x2n-d2;
+% %     c
+% %     d
+%     
+%     p = -0.5;
+% 
+%     in = 0;
+% 
+%             for k=1:length(x2n)
+%                 for l = 1:length(x1n)
+%                     
+%     [my(1) mk(1)] = mex_G00(p,r,u,x1n(l),x2n(k),0);
+%     [my(2) mk(2)] = mex_G00(p,t,v,x1n(l),x2n(k),0);
+%     [my(3) mk(3)] = mex_G00(p,t,u,x1n(l),x2n(k),0);
+%     [my(4) mk(4)] = mex_G00(p,r,v,x1n(l),x2n(k),0);
+% 
+% 
+% %     my_sol = my(1)+my(2)-my(3)-my(4);
+%                     
+%                     in = in + x2w(k) * x1w(l) * (my(1)+my(2)-my(3)-my(4));
+%                     
+%                 end
+%             end
+% %             in
+% 
+% end
\ No newline at end of file
diff --git a/src/test_Fort.m b/src/test_Fort.m
new file mode 100644
index 0000000..ad326d2
--- /dev/null
+++ b/src/test_Fort.m
@@ -0,0 +1,19 @@
+function [ret] = test_Fort(a,b,c,d,r,t,u,v,l1,l2,l3)
+
+
+    quad_sol = surfDoubleQuad(@(x1,x2,y2,y3)1/sqrt((x1-l1).^2+(x2-y2-l2).^2+(y3+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_Fort(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_solveError.m b/src/test_solveError.m
index 9730b85..ae4f5c3 100644
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -1,6 +1,6 @@
 
 %% Laden des Gitters
-load exmpl_2DLShape %Energienorm sollte gegen 8.28466 konvergieren
+load exmpl_3DFichCube %Energienorm sollte gegen 8.28466 konvergieren
 
 %% Rotationen einzelnen von Elementen (zum Testen)
 % coordinates = coordinates(:,[ 1 3 2]);
@@ -16,7 +16,7 @@ end
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
 
-for ref = 1:4
+for ref = 1:3
 
   %% Testen mit Quadratur... und mex_G00
   % A = build_A2(coordinates,elements);