From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Mon, 27 Jun 2011 13:07:18 +0000 (+0000)
Subject: [src] SLPrecangle umbenannt
X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=6a1941f09b8b387cbe36be8911b3663ddecad3db;p=bacc.git

[src] SLPrecangle umbenannt
[src] quadInt in applyInt umbenannt
[doc] beschreibungen in einigen Matlab-Funktionen hinzugefügt

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

diff --git a/src/SLPrecangle.cpp b/src/SLPrecangle.cpp
deleted file mode 100644
index faefd8d..0000000
--- a/src/SLPrecangle.cpp
+++ /dev/null
@@ -1,274 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include <stdlib.h>
-
-#include "SLPrecangle.hpp"
-
-#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 pt = p + 1.5;
-	double sol = 0;
-	if (pt == 0) {
-		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;
-		cout << "warning in G00: no case for p="<< p <<" defined. Possible: [-1.5,-0.5,0.5,...]\n";
-		//mexErrMsgTxt("no case for p defined\n");
-	}
-
-	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);
-	/*
-	 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);
-
-}
-
-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:
-			cout << "p must be either 0, -1/2, -1 or -3/2. (" << p << ")\n";
-			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;
-}
-
-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
deleted file mode 100644
index d801cfe..0000000
--- a/src/SLPrecangle.hpp
+++ /dev/null
@@ -1,34 +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,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(
-		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/build_A.cpp b/src/build_A.cpp
deleted file mode 100644
index 7268af9..0000000
--- a/src/build_A.cpp
+++ /dev/null
@@ -1,289 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-//#include "tbb/parallel_for.h"
-
-
-#include "SLPrecangle.hpp"
-
-#define EPS 0.02
-
-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(Mx4))");
-	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 != 4)
-		mexErrMsgTxt("expected elements (Mx4)");
-	//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 * x0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
-	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 * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
-	double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
-
-	double * y0 = 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 * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
-
-	double tmp;
-
-	//LageInformationen
-	int rx, rxa, rxb, ry, rya, ryb;
-
-	//Ausrechnen
-	for (j = 0; j < em; ++j) {
-		x0[0] = C[(int) E[j] - 1];
-		x0[1] = C[cm + (int) E[j] - 1];
-		x0[2] = C[2 * cm + (int) E[j] - 1];
-
-		x1[0] = C[(int) E[em + j] - 1];
-		x1[1] = C[cm + (int) E[em + j] - 1];
-		x1[2] = C[2 * cm + (int) E[em + j] - 1];
-
-		x2[0] = C[(int) E[2 * em + j] - 1];
-		x2[1] = C[cm + (int) E[2 * em + j] - 1];
-		x2[2] = C[2 * cm + (int) E[2 * em + j] - 1];
-
-		x3[0] = C[(int) E[3 * em + j] - 1];
-		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)
-			xa[i] = x1[i] - x0[i];
-
-		for (i = 0; i < 3; ++i)
-			xb[i] = x3[i] - x0[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]);
-
-		// Lageeigenschaften des Flaechenstuecks
-		//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;
-
-		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) {
-				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];
-				//printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");
-			}
-		}
-
-		//for (i=0;i<3;++i)
-		//	dt[i] = 0;
-
-
-		// printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);
-
-		for (k = 0; k < em; ++k) {
-			y0[0] = C[(int) E[k] - 1];
-			y0[1] = C[cm + (int) E[k] - 1];
-			y0[2] = C[2 * cm + (int) E[k] - 1];
-
-			y1[0] = C[(int) E[em + k] - 1];
-			y1[1] = C[cm + (int) E[em + k] - 1];
-			y1[2] = C[2 * cm + (int) E[em + k] - 1];
-
-			y2[0] = C[(int) E[2 * em + k] - 1];
-			y2[1] = C[cm + (int) E[2 * em + k] - 1];
-			y2[2] = C[2 * cm + (int) E[2 * em + k] - 1];
-
-			y3[0] = C[(int) E[3 * em + k] - 1];
-			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)
-				ya[i] = y1[i] - y0[i];
-
-			for (i = 0; i < 3; ++i)
-				yb[i] = y3[i] - y0[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;
-			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) {
-					for (i = 0; i < 3; ++i)
-						d[i] = dt[i] + y0[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)
-						d[i] = dt[i] + y1[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];
-
-			//printf("(%d,%d)",rx,ry);
-			if (rx == ry) {
-				//printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
-				//printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);
-				//printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);
-				//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
-				if (rxa == rya) {
-					//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]);
-					//	 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,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]);
-					//	 printf("%d%d|%.2f\n",j,k,tmp);
-				}
-
-			} else {
-				//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 * M_PI) * tmp;
-
-		}
-
-	}
-	//Rueckgabe (eventuell zurueck schreiben)
-	//   mxFree(x0);
-	//mxFree(x1);
-	//mxFree(x3);
-	//mxFree(y0);
-	//mxFree(y1);
-	//mxFree(y3);
-	//mxFree(xn);
-	//mxFree(yn);
-	//mxFree(d);
-
-	return;
-}
diff --git a/src/mex_Fort.cpp b/src/mex_Fort.cpp
index 7d4f894..0f4b141 100644
--- a/src/mex_Fort.cpp
+++ b/src/mex_Fort.cpp
@@ -4,7 +4,7 @@
 #include "mex.h"
 #include <stdlib.h>
 
-#include "SLPrecangle.hpp"
+#include "slpRectangle.hpp"
 
 #define my_PI 3.141592653589793
 #define EPS 0.02
diff --git a/src/mex_Fpar.cpp b/src/mex_Fpar.cpp
index ab8ae94..c02cfb6 100644
--- a/src/mex_Fpar.cpp
+++ b/src/mex_Fpar.cpp
@@ -4,7 +4,7 @@
 #include "mex.h"
 #include <stdlib.h>
 
-#include "SLPrecangle.hpp"
+#include "slpRectangle.hpp"
 
 #define my_PI 3.141592653589793
 #define EPS 0.02
diff --git a/src/mex_G00.cpp b/src/mex_G00.cpp
index 0bc77f3..9f9fd98 100644
--- a/src/mex_G00.cpp
+++ b/src/mex_G00.cpp
@@ -4,7 +4,7 @@
 #include "mex.h"
 #include <stdlib.h>
 
-#include "SLPrecangle.hpp"
+#include "slpRectangle.hpp"
 
 #define my_PI 3.141592653589793
 #define EPS 0.02
diff --git a/src/mex_build_A.cpp b/src/mex_build_A.cpp
new file mode 100644
index 0000000..48a30b0
--- /dev/null
+++ b/src/mex_build_A.cpp
@@ -0,0 +1,289 @@
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+//#include "tbb/parallel_for.h"
+
+
+#include "slpRectangle.hpp"
+
+#define EPS 0.02
+
+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(Mx4))");
+	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 != 4)
+		mexErrMsgTxt("expected elements (Mx4)");
+	//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 * x0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	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 * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+
+	double * y0 = 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 * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+
+	double tmp;
+
+	//LageInformationen
+	int rx, rxa, rxb, ry, rya, ryb;
+
+	//Ausrechnen
+	for (j = 0; j < em; ++j) {
+		x0[0] = C[(int) E[j] - 1];
+		x0[1] = C[cm + (int) E[j] - 1];
+		x0[2] = C[2 * cm + (int) E[j] - 1];
+
+		x1[0] = C[(int) E[em + j] - 1];
+		x1[1] = C[cm + (int) E[em + j] - 1];
+		x1[2] = C[2 * cm + (int) E[em + j] - 1];
+
+		x2[0] = C[(int) E[2 * em + j] - 1];
+		x2[1] = C[cm + (int) E[2 * em + j] - 1];
+		x2[2] = C[2 * cm + (int) E[2 * em + j] - 1];
+
+		x3[0] = C[(int) E[3 * em + j] - 1];
+		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)
+			xa[i] = x1[i] - x0[i];
+
+		for (i = 0; i < 3; ++i)
+			xb[i] = x3[i] - x0[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]);
+
+		// Lageeigenschaften des Flaechenstuecks
+		//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;
+
+		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) {
+				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];
+				//printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");
+			}
+		}
+
+		//for (i=0;i<3;++i)
+		//	dt[i] = 0;
+
+
+		// printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);
+
+		for (k = 0; k < em; ++k) {
+			y0[0] = C[(int) E[k] - 1];
+			y0[1] = C[cm + (int) E[k] - 1];
+			y0[2] = C[2 * cm + (int) E[k] - 1];
+
+			y1[0] = C[(int) E[em + k] - 1];
+			y1[1] = C[cm + (int) E[em + k] - 1];
+			y1[2] = C[2 * cm + (int) E[em + k] - 1];
+
+			y2[0] = C[(int) E[2 * em + k] - 1];
+			y2[1] = C[cm + (int) E[2 * em + k] - 1];
+			y2[2] = C[2 * cm + (int) E[2 * em + k] - 1];
+
+			y3[0] = C[(int) E[3 * em + k] - 1];
+			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)
+				ya[i] = y1[i] - y0[i];
+
+			for (i = 0; i < 3; ++i)
+				yb[i] = y3[i] - y0[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;
+			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) {
+					for (i = 0; i < 3; ++i)
+						d[i] = dt[i] + y0[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)
+						d[i] = dt[i] + y1[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];
+
+			//printf("(%d,%d)",rx,ry);
+			if (rx == ry) {
+				//printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
+				//printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);
+				//printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);
+				//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
+				if (rxa == rya) {
+					//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
+							= apply0Int(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 {
+					//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
+							= apply0Int(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);
+				}
+
+			} else {
+				//printf("(%d,%d@%d,%d)", rxa, rxb, rya, ryb);
+				//printf("(%d,%d)", rx, ry);
+				//printf("\n");
+
+				if (rxa == rya) {
+					tmp
+							= apply0Int(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
+							= apply0Int(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
+							= apply0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(ya[rya]), fabs(yb[ryb]), d[rxa],
+									d[rxb], d[rx]);
+				} else {
+					tmp
+							= apply0Int(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 * M_PI) * tmp;
+
+		}
+
+	}
+	//Rueckgabe (eventuell zurueck schreiben)
+	//   mxFree(x0);
+	//mxFree(x1);
+	//mxFree(x3);
+	//mxFree(y0);
+	//mxFree(y1);
+	//mxFree(y3);
+	//mxFree(xn);
+	//mxFree(yn);
+	//mxFree(d);
+
+	return;
+}
diff --git a/src/mex_g0.cpp b/src/mex_g0.cpp
index a93de87..9924813 100644
--- a/src/mex_g0.cpp
+++ b/src/mex_g0.cpp
@@ -4,7 +4,7 @@
 #include "mex.h"
 #include <stdlib.h>
 
-#include "SLPrecangle.hpp"
+#include "slpRectangle.hpp"
 
 #define my_PI 3.141592653589793
 #define EPS 0.02
diff --git a/src/mexit.m b/src/mexit.m
index f552eff..ed3eaad 100644
--- a/src/mexit.m
+++ b/src/mexit.m
@@ -1,8 +1,8 @@
 
-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 mex_g0.cpp slpRectangle.cpp
+mex mex_G00.cpp slpRectangle.cpp
+mex mex_Fpar.cpp slpRectangle.cpp
+mex mex_Fort.cpp slpRectangle.cpp
 
-mex build_A.cpp SLPrecangle.cpp
+mex mex_build_A.cpp slpRectangle.cpp
 
diff --git a/src/quadNorm.m b/src/quadNorm.m
index c663ae4..18b0f85 100644
--- a/src/quadNorm.m
+++ b/src/quadNorm.m
@@ -6,6 +6,7 @@ function n = quadNorm(coordinates, elements,varargin)
 % Diese Funktion Berechnet die Orthogonalen mit Laenge 1 über alle Flächen
 % FLAG:
 % w -> Laenge entspricht Flaecheninhalt
+%
 % P.Schaefer
 
 %% Parameterueberpruefung
diff --git a/src/refineQuad.m b/src/refineQuad.m
index 8502634..aa4c2ae 100644
--- a/src/refineQuad.m
+++ b/src/refineQuad.m
@@ -1,4 +1,18 @@
 function [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type)
+%
+%  [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type)
+%
+% Verfeinert die markierten Elemente mit dem entsprechenden TYP und gibt
+% auch die F2S Beziehungen zurück. type muss von der Länge der Anzahl der
+% Elemente entsprechen und die Einträge können 1,2,3,4 sein.
+%
+% der Typ zu jedem Element entspricht dabei:
+% 1- keine Verfeinerung
+% 2- 4 neue Elemente
+% 3- 2 neue Elemente, übereinander
+% 4- 2 neue Elemente, nebeneinander
+%
+% P. Schaefer
 
 if([1 1] == size(type))
     type = repmat(type, size(elements,1),1);
diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp
new file mode 100644
index 0000000..7e9a7b0
--- /dev/null
+++ b/src/slpRectangle.cpp
@@ -0,0 +1,274 @@
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include <stdlib.h>
+
+#include "slpRectangle.hpp"
+
+#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 pt = p + 1.5;
+	double sol = 0;
+	if (pt == 0) {
+		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;
+		cout << "warning in G00: no case for p="<< p <<" defined. Possible: [-1.5,-0.5,0.5,...]\n";
+		//mexErrMsgTxt("no case for p defined\n");
+	}
+
+	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 applyInt(
+		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 apply0Int(
+		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);
+
+}
+
+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:
+			cout << "p must be either 0, -1/2, -1 or -3/2. (" << p << ")\n";
+			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;
+}
+
+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/slpRectangle.hpp b/src/slpRectangle.hpp
new file mode 100644
index 0000000..26ebd71
--- /dev/null
+++ b/src/slpRectangle.hpp
@@ -0,0 +1,34 @@
+#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);
+
+
+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 applyInt(
+		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 apply0Int(
+		double(*)(double, double, double, double, double, double, double),
+		double, double, double, double, double, double, double);
+
+
+#endif
diff --git a/src/surfDoubleQuad.m b/src/surfDoubleQuad.m
index 86ed43c..4759b7a 100644
--- a/src/surfDoubleQuad.m
+++ b/src/surfDoubleQuad.m
@@ -1,51 +1,12 @@
-% function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
-% %Vierfachintegral ueber f(x1,x2,y1,y2) auf [a b]x[c d]x[r t]x[u v] mit 
-% %der Simpsonregel und w Auswertungspunkten
-%     s = (b-a)/w;
-%     in = 0;
-%     for i = a:s:b-s
-%         in = in + fx2(f,i,c,d,r,t,u,v,w) + ...
-%                 4*fx2(f,(2*i+s)/2,c,d,r,t,u,v,w) + ...
-%                   fx2(f,i+s,c,d,r,t,u,v,w);
-%     end
-%     in = in *s/6;
-% end
-% 
-% function in = fx2(f,x1,c,d,r,t,u,v,w)
-%     s = (d-c)/w;
-%     in = 0;
-%     for i = c:s:d-s
-%         in = in + fy1(f,x1,i,r,t,u,v,w) + ...
-%                 4*fy1(f,x1,(2*i+s)/2,r,t,u,v,w) + ...
-%                   fy1(f,x1,i+s,r,t,u,v,w);
-%     end
-%     in = in *s/6;
-% end
-% 
-% function in = fy1(f,x1,x2,r,t,u,v,w)
-%     s = (t-r)/w;
-%     in = 0;
-%     for i = r:s:t-s
-%         in = in + fy2(f,x1,x2,i,u,v,w) + ...
-%                 4*fy2(f,x1,x2,(2*i+s)/2,u,v,w) + ...
-%                   fy2(f,x1,x2,i+s,u,v,w);
-%     end
-%     in = in *s/6;
-% end
-% 
-% function in = fy2(f,x1,x2,y1,u,v,w)
-%     s = (v-u)/w;
-%     in = 0;
-%     for i = u:s:v-s
-%         in = in + f(x1,x2,y1,i) + ...
-%                 4*f(x1,x2,y1,(2*i+s)/2) + ...
-%                   f(x1,x2,y1,i+s);
-%     end
-%     in = in *s/6;
-% end
-
-
 function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
+%
+% in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
+%
+% Berechnet das Vierfachintegral über die Funktion f(x1,x2,y1,y2) mit den 
+% Grenzen [a b]x[c d]x[r t]x[u v]. Dazu wird eine Gaussquadratur über w 
+% Punkte verwendet.
+%
+% P. Schaefer
 
     [x1n x1w] = gauss(w,a,b);
     [x2n x2w] = gauss(w,c,d);
@@ -58,47 +19,12 @@ function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
             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));
+                    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;
-% 
-%     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
+end
\ No newline at end of file
diff --git a/src/surfQuad.m b/src/surfQuad.m
index 3111d7d..36b2289 100644
--- a/src/surfQuad.m
+++ b/src/surfQuad.m
@@ -1,25 +1,12 @@
-% function in = surfQuad(f,a,b,c,d,v)
-% %Doppelintegral ueber f(x,y) auf [a b]x[c d] mit der Simpsonregel und v
-% %Auswertungspunkten
-%     s = (b-a)/v;
-%     in = 0;
-%     for i = a:s:b-s
-%         in = in + fy(f,i,c,d,v) + 4*fy(f,(2*i+s)/2,c,d,v) + fy(f,i+s,c,d,v);
-%     end
-%     in = in *s/6;
-% end
-% 
-% function in = fy(f,x,c,d,v)
-%     s = (d-c)/v;
-%     in = 0;
-%     for i = c:s:d-s
-%         in = in + f(x,i) + 4*f(x,(2*i+s)/2) + f(x,i+s);
-%     end
-%     in = in *s/6;
-% end
-
-
 function in = surfQuad(f,a,b,c,d,v)
+%
+% in = surfQuad(f,a,b,c,d,v)
+%
+% Berechnet das Doppelintegral über die Funktion f(x,y) mit den Grenzen
+% [a b]x[c d]. Dazu wird eine Gaussquadratur über v Punkte verwendet.
+%
+% P. Schaefer
+
     [xn xw] = gauss(v,a,b);
     [yn yw] = gauss(v,c,d);
     
diff --git a/src/test_solve.m b/src/test_solve.m
index 184dc97..c6245c7 100644
--- a/src/test_solve.m
+++ b/src/test_solve.m
@@ -2,7 +2,7 @@
 N = size(elements,1);
 
 tic
-A = build_A(coordinates,elements);
+A = mex_build_A(coordinates,elements);
 t = toc;
 b = sqrt(sum(quadNorm(coordinates,elements,'w').^2,2));
 
diff --git a/src/test_solveError.m b/src/test_solveError.m
index fa16e15..b9a94a5 100644
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -16,7 +16,7 @@ while size(elements,1)<maxSize
   ref = ref+1;
 
   [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
-  A_fine = build_A(coordinates_fine,elements_fine);
+  A_fine = mex_build_A(coordinates_fine,elements_fine);
   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
   x_fine = A_fine\b;
 %   xe_fine = x_fine'*A_fine*x_fine;
@@ -41,7 +41,7 @@ while size(elements,1)<maxSize
   ref = ref+1;
 
   [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
-  A_fine = build_A(coordinates_fine,elements_fine);
+  A_fine = mex_build_A(coordinates_fine,elements_fine);
   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
   x_fine = A_fine\b;
 %   xe_fine = x_fine'*A_fine*x_fine;
diff --git a/src/test_solveS.m b/src/test_solveS.m
index 1ea205e..96d646e 100644
--- a/src/test_solveS.m
+++ b/src/test_solveS.m
@@ -18,7 +18,7 @@ A2 = build_A2(coordinates,elements);
 t(1) = toc;
 disp ' '
 tic
-A2 = build_A(coordinates,elements);
+A2 = mex_build_A(coordinates,elements);
 t(2) = toc;
 b = sqrt(sum(quadNorm(coordinates,elements,'w').^2,2));