From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Thu, 20 Oct 2011 19:06:32 +0000 (+0000)
Subject: [src] neue Funktionen zum Automatisieren der Anisotropen und Isotropen Berechnung A_*
X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=afcb2750a180a998a2ae0862656678a2958341e4;p=bacc.git

[src] neue Funktionen zum Automatisieren der Anisotropen und Isotropen Berechnung A_*
[src] nicht verwendete Skripte entfernt
[src] slpRect erweitert
[src] mex_build_A ist nun von Parametern abhängig die bestimmen wie Berechnet werden soll

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

diff --git a/src/A_loadMesh.m b/src/A_loadMesh.m
new file mode 100644
index 0000000..703edb4
--- /dev/null
+++ b/src/A_loadMesh.m
@@ -0,0 +1,11 @@
+function A_loadMesh(coordinates,elements)
+% Laed ein Mesh
+%
+% A_loadMesh(coordinates,elements)
+
+  global G_C;
+  global G_E;
+  G_C = coordinates;
+  G_E = elements;
+end
+
diff --git a/src/A_setParam.m b/src/A_setParam.m
new file mode 100644
index 0000000..30cc2a2
--- /dev/null
+++ b/src/A_setParam.m
@@ -0,0 +1,8 @@
+function A_setParam(paras)
+%setzt die Parameter fuer die Durchlaeufe
+
+global G_para;
+G_para = paras;
+
+end
+
diff --git a/src/A_stepAnIso.m b/src/A_stepAnIso.m
new file mode 100644
index 0000000..5e4d976
--- /dev/null
+++ b/src/A_stepAnIso.m
@@ -0,0 +1,30 @@
+function dataAniso = A_stepAnIso(mu,type)
+
+global G_E;
+global G_C;
+
+if (isempty(G_E) || isempty(G_C))
+  disp 'Error: Please use A_loadMesh first'
+  return
+end
+
+  [coordinates_fine,elements_fine, f2s]=refineQuad(G_C,G_E,2*ones(1,size(G_E,1)));
+  
+  A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type);
+  
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+  
+  x_fine = A_fine\b;
+  
+  xe_fine = x_fine'*A_fine*x_fine;
+
+  ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
+  
+  marked = mark(x_fine(f2s)',ind,0.5,0.5);
+  
+  dataAniso = [size(G_E,1) sqrt(sum(ind)) xe_fine];
+  [G_C, G_E] = refineQuad(G_C,G_E,marked);
+
+
+end
+
diff --git a/src/A_stepIso.m b/src/A_stepIso.m
new file mode 100644
index 0000000..264f245
--- /dev/null
+++ b/src/A_stepIso.m
@@ -0,0 +1,28 @@
+function dataIso = A_stepIso(mu,type)
+
+global G_E;
+global G_C;
+
+if (isempty(G_E) || isempty(G_C))
+  disp 'Error: Please use A_loadMesh first'
+  return
+end
+
+  [coordinates_fine,elements_fine, f2s]=refineQuad(G_C,G_E,2*ones(1,size(G_E,1)));
+  
+  A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type);
+  
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+  
+  x_fine = A_fine\b;
+  
+  xe_fine = x_fine'*A_fine*x_fine;
+
+  ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
+  
+  dataIso = [size(G_E,1) sqrt(sum(ind)) xe_fine];
+  [G_C, G_E] = refineQuad(G_C,G_E,2*ones(1,size(G_E,1)));
+
+
+end
+
diff --git a/src/build_A2.m b/src/build_A2.m
deleted file mode 100644
index 9550b53..0000000
--- a/src/build_A2.m
+++ /dev/null
@@ -1,35 +0,0 @@
-function A = build_A2(coordinates,elements)
-    N = size(elements,1);
-
-    A1 = zeros(N);
-
-    % untere schranke s t obere schranke k l
-    intF = @(f,a,b,c,d,r,t,u,v)...
-        f(b,d,t,v)-f(b,d,t,u)-f(b,d,r,v)+f(b,d,r,u)...
-        -f(b,c,t,v)+f(b,c,t,u)+f(b,c,r,v)-f(b,c,r,u)...
-        -f(a,d,t,v)+f(a,d,t,u)+f(a,d,r,v)-f(a,d,r,u)...
-        +f(a,c,t,v)-f(a,c,t,u)-f(a,c,r,v)+f(a,c,r,u);
-    
-    for j = 1:N
-        for k = 1:N
-            ej = coordinates(elements(j,:)',:);
-            ek = coordinates(elements(k,:)',:);
-            
-            d = ek(1,:) - ej(1,:);
-            
-            ej = ej - repmat(ej(1,:),3,1);
-            ek = ek - repmat(ek(1,:),3,1);
-
-%             d = ones(1,3);
-
-%             A1(j,k) = surfDoubleQuad(@(x1,x2,y1,y2) 1/sqrt((x1-y1-d(1)).^2+(x2-y2-d(2)).^2+d(3).^2)...
-%                 ,ej(1,1),ej(2,1),ej(1,2), ej(3,2),ek(1,1), ek(2,1),ek(1,2), ek(3,2),4);
-
-            A1(j,k) = intF(@(x1,x2,y1,y2) mex_Fpar(x1,x2,y1,y2,d(1),d(2),d(3))...
-                ,ej(1,1),ej(2,1),ej(1,2), ej(3,2),ek(1,1), ek(2,1),ek(1,2), ek(3,2));
-        end
-    end
-    
-    A = A1/(4*pi);
-
-end
\ No newline at end of file
diff --git a/src/mex_Fort.cpp b/src/mex_Fort.cpp
deleted file mode 100644
index 4741cb0..0000000
--- a/src/mex_Fort.cpp
+++ /dev/null
@@ -1,30 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.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);
-
-	*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]));
-
-}
-
-
-
diff --git a/src/mex_Fpar.cpp b/src/mex_Fpar.cpp
deleted file mode 100644
index 369f683..0000000
--- a/src/mex_Fpar.cpp
+++ /dev/null
@@ -1,30 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.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,y1,y2,d1,d2,d3)");
-	if (nlhs > 2)
-		mexErrMsgTxt("has maximal two output arguments");
-
-
-	plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
-	*mxGetPr(plhs[0]) = F_par(*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_G00.cpp b/src/mex_G00.cpp
deleted file mode 100644
index 1c946f8..0000000
--- a/src/mex_G00.cpp
+++ /dev/null
@@ -1,32 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.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 != 6)
-		mexErrMsgTxt("expected (p, y1, y2, x1, x2, l)");
-	if (nlhs > 2)
-		mexErrMsgTxt("has maximal two output arguments");
-
-
-	plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
-	plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
-	*mxGetPr(plhs[0]) = G_AY1Y2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
-	*mxGetPr(plhs[1]) = G_QY1Y2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
-
-}
-
-
-
diff --git a/src/mex_build_A.cpp b/src/mex_build_A.cpp
deleted file mode 100644
index 5cd019f..0000000
--- a/src/mex_build_A.cpp
+++ /dev/null
@@ -1,296 +0,0 @@
-/*
- * voll Analytisch
- * ohne Optimierungen
- *
- *
- */
-
-
-#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
-							= calcParIntA( 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
-							= calcParIntA( 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
-							= calcOrtIntA( fabs(xb[rxb]), fabs(xa[rxa]),
-									fabs(ya[rya]), fabs(yb[ryb]), d[rxb],
-									d[rxa], d[rx]);
-				} else if (rxa == ryb) {
-					tmp
-							= calcOrtIntA( fabs(xb[rxb]), fabs(xa[rxa]),
-									fabs(yb[ryb]), fabs(ya[rya]), d[rxb],
-									d[rxa], d[rx]);
-				} else if (rxb == rya) {
-					tmp
-							= calcOrtIntA( fabs(xa[rxa]), fabs(xb[rxb]),
-									fabs(ya[rya]), fabs(yb[ryb]), d[rxa],
-									d[rxb], d[rx]);
-				} else {
-					tmp
-							= calcOrtIntA( fabs(xa[rxa]), fabs(xb[rxb]),
-									fabs(yb[ryb]), fabs(ya[rya]), d[rxa],
-									d[rxb], d[rx]);
-				}
-
-			}
-			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_build_AU.cpp b/src/mex_build_AU.cpp
new file mode 100644
index 0000000..9519a0d
--- /dev/null
+++ b/src/mex_build_AU.cpp
@@ -0,0 +1,276 @@
+/*
+ * Art der Berechnung durch Parameter bestimmt...
+ * Analytisch/Semianalytisch/...
+ *
+ *
+ */
+
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+#define M_EPS 10e-8
+
+//#include "tbb/parallel_for.h"
+
+
+#include "slpRectangle.hpp"
+
+using namespace std;
+
+int dimOfVec(double* vec) {
+	if (fabs(vec[2]) > M_EPS)
+		return 2;
+	else if (fabs(vec[1]) > M_EPS)
+		return 1;
+	else if (fabs(vec[0]) > M_EPS)
+		return 0;
+	else {
+		cerr << "dimOfVec : (" << vec[0] << vec[1] << vec[2]
+				<< ") alle Werte 0 \n";
+		return -1;
+	}
+
+}
+
+inline int dimOfThird(int a, int b) {
+	return ((-(a + b) % 3) + 3) % 3;
+}
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+
+	int i, j, k; //Schleifenindizes
+	double tmp; //Zwischenspeicherung der Einzelnen Werte
+
+	//Sicherheitsabfragen zu Datengroessen
+	if ((nrhs != 4))
+		mexErrMsgTxt(
+				"expected (coordinates(Nx3),elements(Mx4),mu(double),type(int))");
+	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)");
+	cout << "  Coordinaten:" << cm << endl;
+
+	int em = mxGetM(prhs[1]);
+	int en = mxGetN(prhs[1]);
+	if (en != 4)
+		mexErrMsgTxt("expected elements (Mx4)");
+	cout << "  Elemente:" << em << endl;
+
+	//Auslesen der Parameter
+
+	plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);
+	double * A = mxGetPr(plhs[0]);
+	double * C = mxGetPr(prhs[0]);
+	double * E = mxGetPr(prhs[1]);
+
+	double mu = *(mxGetPr(prhs[2]));
+	int type = *(mxGetPr(prhs[3]));
+
+	//Initialisieren der Hilfsvektoren
+
+	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 * 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 * 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));
+
+	// Welche Funktion soll verwendet werden
+	double(*ctypeP)(double, double, double, double, double, double, double,
+			double);
+	double(*ctypeO)(double, double, double, double, double, double, double,
+			double);
+
+	//Art der Berechnung bestimmen
+	cout << "  Type:" << type << endl;
+	switch (type) {
+	case 0:
+		ctypeP = calcParIntO0;
+		ctypeO = calcOrtIntO0;
+	case 1:
+		ctypeP = calcParIntO1;
+		ctypeO = calcOrtIntO1;
+	case 2:
+		ctypeP = calcParIntO2;
+		ctypeO = calcOrtIntO2;
+	case 3:
+		ctypeP = calcParIntO3;
+		ctypeO = calcOrtIntO3;
+
+	}
+
+	//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];
+
+		// Lageeigenschaften des Flaechenstuecks
+
+		rxa = dimOfVec(xa);
+		rxb = dimOfVec(xb);
+		rx = dimOfThird(rxa, rxb);
+
+		//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];
+			}
+		}
+
+		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];
+
+			// Lageeigenschaften des Flaechenstuecks
+
+			rya = dimOfVec(ya);
+			ryb = dimOfVec(yb);
+			ry = dimOfThird(rya, ryb);
+
+			//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];
+				}
+			}
+
+			if (rx == ry) {
+				if (rxa == rya) {
+					tmp
+							= ctypeP(fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(ya[rxa]), fabs(yb[rxb]), d[rxa],
+									d[rxb], d[rx], mu);
+
+				} else {
+					tmp
+							= ctypeP(fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(yb[rxa]), fabs(ya[rxb]), d[rxa],
+									d[rxb], d[rx], mu);
+				}
+
+			} else {
+
+				if (rxa == rya) {
+					tmp
+							= ctypeO(fabs(xb[rxb]), fabs(xa[rxa]),
+									fabs(ya[rya]), fabs(yb[ryb]), d[rxb],
+									d[rxa], d[rx], mu);
+				} else if (rxa == ryb) {
+					tmp
+							= ctypeO(fabs(xb[rxb]), fabs(xa[rxa]),
+									fabs(yb[ryb]), fabs(ya[rya]), d[rxb],
+									d[rxa], d[rx], mu);
+				} else if (rxb == rya) {
+					tmp
+							= ctypeO(fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(ya[rya]), fabs(yb[ryb]), d[rxa],
+									d[rxb], d[rx], mu);
+				} else {
+					tmp
+							= ctypeO(fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(yb[ryb]), fabs(ya[rya]), d[rxa],
+									d[rxb], d[rx], mu);
+				}
+
+			}
+			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(d);
+	//mxFree(dt);
+
+	return;
+}
diff --git a/src/mex_build_As1.cpp b/src/mex_build_As1.cpp
deleted file mode 100644
index cc1b7e5..0000000
--- a/src/mex_build_As1.cpp
+++ /dev/null
@@ -1,287 +0,0 @@
-/*
- * voll Analytisch
- * Element mit groesserer Diagonale nach aussen gedreht
- *
- *
- */
-
-
-#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;
-
-inline int dimOfVec(double * vec)
-{
-	if (vec[2] != 0)
-		return 2;
-	else if (vec[1] != 0)
-		return 1;
-	else
-		return 0;
-}
-
-
-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,tmp2;
-
-	int itmp;
-	double dtmp;
-	double * vtmp;
-
-	long double xda,xdb,yda,ydb;
-
-	//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])
-
-		rx = dimOfVec(xn);
-
-		rxa = dimOfVec(xa);
-
-		rxb = dimOfVec(xb);
-
-
-		//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");
-			}
-		}
-
-		xda = xa[rxa];
-		xdb = xb[rxb];
-
-		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]);
-
-			ry = dimOfVec(yn);
-
-			rya = dimOfVec(ya);
-
-			ryb = dimOfVec(yb);
-
-			//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");
-				}
-			}
-
-			yda = ya[rya];
-			ydb = yb[ryb];
-
-			//printf("(%d,%d)",rx,ry);
-			if (rx == ry) {
-				//printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
-				//printf("[%.2f,%.2f@%.2f,%.2f]",xda,xdb,yda,ydb);
-				//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
-							= calcParIntO1( fabs(xda), fabs(xdb),
-									fabs(yda), fabs(ydb), d[rxa],
-									d[rxb], d[rx]);
-
-				} 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
-							= calcParIntO1( fabs(xda), fabs(xdb),
-									fabs(ydb), fabs(yda), 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
-							= calcOrtIntO1( fabs(xdb), fabs(xda),
-									fabs(yda), fabs(ydb), d[rxb],
-									d[rxa], d[rx]);
-				} else if (rxa == ryb) {
-					tmp
-							= calcOrtIntO1( fabs(xdb), fabs(xda),
-									fabs(ydb), fabs(yda), d[rxb],
-									d[rxa], d[rx]);
-				} else if (rxb == rya) {
-					tmp
-							= calcOrtIntO1( fabs(xda), fabs(xdb),
-									fabs(yda), fabs(ydb), d[rxa],
-									d[rxb], d[rx]);
-				} else {
-					tmp
-							= calcOrtIntO1( fabs(xda), fabs(xdb),
-									fabs(ydb), fabs(yda), d[rxa],
-									d[rxb], d[rx]);
-				}
-
-			}
-			A[(k * em) + j] = 1. / (4 * M_PI) * tmp;
-
-		}
-
-	}
-	//printf("\n");
-	//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_build_As2.cpp b/src/mex_build_As2.cpp
deleted file mode 100644
index 2366177..0000000
--- a/src/mex_build_As2.cpp
+++ /dev/null
@@ -1,298 +0,0 @@
-/*
- * voll Analytisch
- * Element mit groesserer Diagonale nach aussen gedreht
- * Quadratur ueber E_j bei dist(E_j,E_k) >= eta*dia(E_j) | dia(E_j) <= dia(E_k)
- *
- */
-
-
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-//#include "tbb/parallel_for.h"
-
-
-#include "slpRectangle.hpp"
-
-//Quadraturverwendung!!!
-#define eta 0.8
-
-using namespace std;
-
-inline int dimOfVec(double * vec)
-{
-
-	if (vec[2] != 0)
-		return 2;
-	else if (vec[1] != 0)
-		return 1;
-	else if (vec[0] != 0)
-		return 0;
-	else{
-		cerr << "dimOfVec : alle Werte 0";
-		return -1;
-	}
-}
-
-inline int dimOfThird(int a, int b)
-{
-	return mod(-(a+b),3);
-}
-
-
-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,tmp2;
-
-	int itmp;
-	double dtmp;
-	double * vtmp;
-
-	long double xda,xdb,yda,ydb;
-
-	//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])
-
-		rx = dimOfVec(xn);
-
-		rxa = dimOfVec(xa);
-
-		rxb = dimOfVec(xb);
-
-
-		//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");
-			}
-		}
-
-		xda = xa[rxa];
-		xdb = xb[rxb];
-
-		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]);
-
-			ry = dimOfVec(yn);
-
-			rya = dimOfVec(ya);
-
-			ryb = dimOfVec(yb);
-
-			//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");
-				}
-			}
-
-			yda = ya[rya];
-			ydb = yb[ryb];
-
-			//printf("(%d,%d)",rx,ry);
-			if (rx == ry) {
-				//printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
-				//printf("[%.2f,%.2f@%.2f,%.2f]",xda,xdb,yda,ydb);
-				//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
-							= calcParIntO2( fabs(xda), fabs(xdb),
-									fabs(yda), fabs(ydb), d[rxa],
-									d[rxb], d[rx],eta);
-
-				} 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
-							= calcParIntO2( fabs(xda), fabs(xdb),
-									fabs(ydb), fabs(yda), d[rxa],
-									d[rxb], d[rx],eta);
-					//	 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
-							= calcOrtIntO2( fabs(xdb), fabs(xda),
-									fabs(yda), fabs(ydb), d[rxb],
-									d[rxa], d[rx],eta);
-				} else if (rxa == ryb) {
-					tmp
-							= calcOrtIntO2( fabs(xdb), fabs(xda),
-									fabs(ydb), fabs(yda), d[rxb],
-									d[rxa], d[rx],eta);
-				} else if (rxb == rya) {
-					tmp
-							= calcOrtIntO2( fabs(xda), fabs(xdb),
-									fabs(yda), fabs(ydb), d[rxa],
-									d[rxb], d[rx],eta);
-				} else {
-					tmp
-							= calcOrtIntO2( fabs(xda), fabs(xdb),
-									fabs(ydb), fabs(yda), d[rxa],
-									d[rxb], d[rx],eta);
-				}
-
-			}
-			A[(k * em) + j] = 1. / (4 * M_PI) * tmp;
-
-		}
-
-	}
-	//printf("\n");
-	//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_build_As3.cpp b/src/mex_build_As3.cpp
deleted file mode 100644
index 6a6da3a..0000000
--- a/src/mex_build_As3.cpp
+++ /dev/null
@@ -1,288 +0,0 @@
-/*
- * voll Analytisch
- * Element mit groesserer Diagonale nach aussen gedreht
- * Quadratur ueber E_j bei dist(E_j,E_k) >= eta*dia(E_j) | dia(E_j) <= dia(E_k)
- *
- */
-
-
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-//#include "tbb/parallel_for.h"
-
-
-#include "slpRectangle.hpp"
-
-//Quadraturverwendung!!!
-#define eta 0.8
-
-using namespace std;
-
-inline int dimOfVec(double * vec)
-{
-	if (vec[2] != 0)
-		return 2;
-	else if (vec[1] != 0)
-		return 1;
-	else
-		return 0;
-}
-
-
-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,tmp2;
-
-	int itmp;
-	double dtmp;
-	double * vtmp;
-
-	long double xda,xdb,yda,ydb;
-
-	//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])
-
-		rx = dimOfVec(xn);
-
-		rxa = dimOfVec(xa);
-
-		rxb = dimOfVec(xb);
-
-
-		//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");
-			}
-		}
-
-		xda = xa[rxa];
-		xdb = xb[rxb];
-
-		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]);
-
-			ry = dimOfVec(yn);
-
-			rya = dimOfVec(ya);
-
-			ryb = dimOfVec(yb);
-
-			//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");
-				}
-			}
-
-			yda = ya[rya];
-			ydb = yb[ryb];
-
-			//printf("(%d,%d)",rx,ry);
-			if (rx == ry) {
-				//printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
-				//printf("[%.2f,%.2f@%.2f,%.2f]",xda,xdb,yda,ydb);
-				//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
-							= calcParIntO3( fabs(xda), fabs(xdb),
-									fabs(yda), fabs(ydb), d[rxa],
-									d[rxb], d[rx],eta);
-
-				} 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
-							= calcParIntO3( fabs(xda), fabs(xdb),
-									fabs(ydb), fabs(yda), d[rxa],
-									d[rxb], d[rx],eta);
-					//	 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
-							= calcOrtIntO3( fabs(xdb), fabs(xda),
-									fabs(yda), fabs(ydb), d[rxb],
-									d[rxa], d[rx],eta);
-				} else if (rxa == ryb) {
-					tmp
-							= calcOrtIntO3( fabs(xdb), fabs(xda),
-									fabs(ydb), fabs(yda), d[rxb],
-									d[rxa], d[rx],eta);
-				} else if (rxb == rya) {
-					tmp
-							= calcOrtIntO3( fabs(xda), fabs(xdb),
-									fabs(yda), fabs(ydb), d[rxa],
-									d[rxb], d[rx],eta);
-				} else {
-					tmp
-							= calcOrtIntO3( fabs(xda), fabs(xdb),
-									fabs(ydb), fabs(yda), d[rxa],
-									d[rxb], d[rx],eta);
-				}
-
-			}
-			A[(k * em) + j] = 1. / (4 * M_PI) * tmp;
-
-		}
-
-	}
-	//printf("\n");
-	//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_calcOrtInt.cpp b/src/mex_calcOrtInt.cpp
deleted file mode 100755
index f4b35ae..0000000
--- a/src/mex_calcOrtInt.cpp
+++ /dev/null
@@ -1,31 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.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,y1,y2,d1,d2,d3)");
-	if (nlhs > 2)
-		mexErrMsgTxt("has maximal two output arguments");
-
-	plhs[0] = mxCreateDoubleMatrix(1, 3, mxREAL);
-
-	*mxGetPr(plhs[0]) = calcortIntA(*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) = calcParIntQX1X2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-//	*(mxGetPr(plhs[0])+2) = calcParIntQY1X1(*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_calcParInt.cpp b/src/mex_calcParInt.cpp
deleted file mode 100755
index 08b26a9..0000000
--- a/src/mex_calcParInt.cpp
+++ /dev/null
@@ -1,32 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.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,y1,y2,d1,d2,d3)");
-	if (nlhs > 2)
-		mexErrMsgTxt("has maximal two output arguments");
-
-
-	plhs[0] = mxCreateDoubleMatrix(1, 3, mxREAL);
-
-	*mxGetPr(plhs[0]) = calcParIntA(*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) = calcParIntQX1X2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-	*(mxGetPr(plhs[0])+2) = calcParIntQY1X1(*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_g0.cpp b/src/mex_g0.cpp
deleted file mode 100644
index 4353069..0000000
--- a/src/mex_g0.cpp
+++ /dev/null
@@ -1,32 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.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 != 4)
-		mexErrMsgTxt("expected (p, y, x, l)");
-	if (nlhs > 2)
-		mexErrMsgTxt("has maximal two output arguments");
-
-
-	plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
-	plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
-	*mxGetPr(plhs[0]) = g_AY(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
-	*mxGetPr(plhs[0]) = g_QY(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
-
-}
-
-
-
diff --git a/src/mexit.m b/src/mexit.m
deleted file mode 100644
index f868508..0000000
--- a/src/mexit.m
+++ /dev/null
@@ -1,10 +0,0 @@
-
-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 mex_build_A.cpp slpRectangle.cpp
-mex mex_build_As1.cpp slpRectangle.cpp
-mex mex_build_As2.cpp slpRectangle.cpp
-
diff --git a/src/save1.mat b/src/save1.mat
deleted file mode 100644
index bc8e6ba..0000000
Binary files a/src/save1.mat and /dev/null differ
diff --git a/src/save10.mat b/src/save10.mat
deleted file mode 100644
index ebe7708..0000000
Binary files a/src/save10.mat and /dev/null differ
diff --git a/src/save11.mat b/src/save11.mat
deleted file mode 100644
index 9ba30c9..0000000
Binary files a/src/save11.mat and /dev/null differ
diff --git a/src/save12.mat b/src/save12.mat
deleted file mode 100644
index 5806d47..0000000
Binary files a/src/save12.mat and /dev/null differ
diff --git a/src/save13.mat b/src/save13.mat
deleted file mode 100644
index 99aadeb..0000000
Binary files a/src/save13.mat and /dev/null differ
diff --git a/src/save2.mat b/src/save2.mat
deleted file mode 100644
index d560b6b..0000000
Binary files a/src/save2.mat and /dev/null differ
diff --git a/src/save3.mat b/src/save3.mat
deleted file mode 100644
index 193c06e..0000000
Binary files a/src/save3.mat and /dev/null differ
diff --git a/src/save4.mat b/src/save4.mat
deleted file mode 100644
index 290ccd8..0000000
Binary files a/src/save4.mat and /dev/null differ
diff --git a/src/save5.mat b/src/save5.mat
deleted file mode 100644
index dde8537..0000000
Binary files a/src/save5.mat and /dev/null differ
diff --git a/src/save6.mat b/src/save6.mat
deleted file mode 100644
index 2fa67a3..0000000
Binary files a/src/save6.mat and /dev/null differ
diff --git a/src/save7.mat b/src/save7.mat
deleted file mode 100644
index bfabafd..0000000
Binary files a/src/save7.mat and /dev/null differ
diff --git a/src/save8.mat b/src/save8.mat
deleted file mode 100644
index 32661e3..0000000
Binary files a/src/save8.mat and /dev/null differ
diff --git a/src/save9.mat b/src/save9.mat
deleted file mode 100644
index 47df2c7..0000000
Binary files a/src/save9.mat and /dev/null differ
diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp
index 8fad411..6a29495 100644
--- a/src/slpRectangle.cpp
+++ b/src/slpRectangle.cpp
@@ -7,14 +7,10 @@
 
 #include "slpRectangle.hpp"
 
-#define EPS 10e5
-
 using namespace std;
 
-double gauss_w5[] = {0.1185,0.2393,0.2844,0.2393,0.1185};
-double gauss_n5[] = {0.0469,0.2308,0.5000,0.7692,0.9531};
-
-
+double gauss_w5[] = { 0.1185, 0.2393, 0.2844, 0.2393, 0.1185 };
+double gauss_n5[] = { 0.0469, 0.2308, 0.5000, 0.7692, 0.9531 };
 
 int sign(double);
 //double arsinh(double);
@@ -22,6 +18,7 @@ int sign(double);
 
 /* ============================== */
 
+
 int inline sign(double x) {
 	return x > 0 ? 1 : (x < 0 ? -1 : 0);
 }
@@ -31,16 +28,17 @@ int inline sign(double x) {
  }
  */
 
-double inline max(double x, double y){
-	return x<y?y:x;
+double inline max(double x, double y) {
+	return x < y ? y : x;
 }
 
 double g_QY(double p, double y, double x, double l) {
 
 	double sol = 0;
 
-	for(int i = 0;i<5;++i){
-		sol += pow((x-gauss_n5[i]*y)*(x-gauss_n5[i]*y)+l*l,p)*gauss_w5[i]*y;
+	for (int i = 0; i < 5; ++i) {
+		sol += pow((x - gauss_n5[i] * y) * (x - gauss_n5[i] * y) + l * l, p)
+				* gauss_w5[i] * y;
 	}
 
 	return sol;
@@ -81,26 +79,27 @@ double g_AY(double p, double y, double x, double l) {
 double G_QY1Y2(double p, double y1, double y2, double x1, double x2, double l) {
 
 	double sol = 0;
-	for(int i = 0;i<5;++i){
-		for(int j=0;j<5;++j){
-			sol += pow((x1-y1*gauss_n5[i])*(x1-y1*gauss_n5[i])
-					+(x2-y2*gauss_n5[j])*(x2-y2*gauss_n5[j])+l*l,p)*y1*gauss_w5[i]*y2*gauss_w5[j];
+	for (int i = 0; i < 5; ++i) {
+		for (int j = 0; j < 5; ++j) {
+			sol += pow(
+					(x1 - y1 * gauss_n5[i]) * (x1 - y1 * gauss_n5[i]) + (x2
+							- y2 * gauss_n5[j]) * (x2 - y2 * gauss_n5[j]) + l
+							* l, p) * y1 * gauss_w5[i] * y2 * gauss_w5[j];
 		}
 	}
 	return sol;
 }
 
-
 double G_AY2X2(double y1, double y2, double x1, double x2, double d3) {
 
 	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);
 
-	double hlp = ((x1 - y1 ) * (x1 - y1 ) + (x2 - y2) * (x2 - y2) + d3 * d3);
+	double hlp = ((x1 - y1) * (x1 - y1) + (x2 - y2) * (x2 - y2) + d3 * d3);
 
 	double sol = sqrt(hlp);
 
 	if ((x2 - y2) != 0)
-		sol += (x2 - y2) * log(sqrt(hlp)-(x2 - y2));
+		sol += (x2 - y2) * log(sqrt(hlp) - (x2 - y2));
 
 	return sol;
 }
@@ -132,13 +131,15 @@ double G_AY1Y2(double p, double y1, double y2, double x1, double x2, double l) {
 		sol /= 2 * p + 2;
 	} else {
 		sol = NAN;
-		cout << "warning in G_AY1Y2: no case for p="<< p <<" defined. Possible: [-1.5,-0.5,0.5,...]\n";
+		cout << "warning in G_AY1Y2: no case for p=" << p
+				<< " defined. Possible: [-1.5,-0.5,0.5,...]\n";
 	}
 
 	return sol;
 }
 
-double Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2, double l) {
+double Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2,
+		double l) {
 	double sol = 0;
 
 	//sol += y2*G_AY1Y2(y1,x2,)
@@ -169,8 +170,6 @@ 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) {
 
@@ -197,19 +196,19 @@ double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,
 }
 
 /*double applyInt4(
-		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) {
+ 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(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 apply0Int4(
 		double(*f)(double, double, double, double, double, double, double),
@@ -229,15 +228,14 @@ double apply0Int2(
 		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, 0, t, v, d1, d2, d3)
-			-f(0, d, t, v, d1, d2, d3)
-			+f(0, 0, t, v, d1, d2, d3);
+	return f(b, d, t, v, d1, d2, d3) - f(b, 0, t, v, d1, d2, d3) - f(0, d, t,
+			v, d1, d2, d3) + f(0, 0, t, v, d1, d2, d3);
 
 }
 
 // Berechnet das eigentliche Integral für parallele Elemente voll Analytisch
-double calcParIntA(double b, double d, double t, double v, double d1, double d2, double d3) {
+double calcParIntA(double b, double d, double t, double v, double d1,
+		double d2, double d3) {
 	return apply0Int4(F_par, b, d, t, v, d1, d2, d3);
 
 }
@@ -271,10 +269,16 @@ double calcParIntQY1X1(double b, double d, double t, double v, double d1,
 
 	for (int i = 0; i < 5; ++i) {
 		for (int j = 0; j < 5; ++j) {
-			sol += G_AY2X2( t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i],d , d3) * t * b * gauss_w5[i] * gauss_w5[j];
-			sol -= G_AY2X2( t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], 0, d3) * t * b * gauss_w5[i] * gauss_w5[j];
-			sol -= G_AY2X2( t * gauss_n5[j] + d1, d2, b * gauss_n5[i], d, d3) * t * b * gauss_w5[i] * gauss_w5[j];
-			sol += G_AY2X2( t * gauss_n5[j] + d1, d2, b * gauss_n5[i], 0, d3) * t * b * gauss_w5[i] * gauss_w5[j];
+			sol
+					+= G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i],
+							d, d3) * t * b * gauss_w5[i] * gauss_w5[j];
+			sol
+					-= G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i],
+							0, d3) * t * b * gauss_w5[i] * gauss_w5[j];
+			sol -= G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], d, d3)
+					* t * b * gauss_w5[i] * gauss_w5[j];
+			sol += G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], 0, d3)
+					* t * b * gauss_w5[i] * gauss_w5[j];
 		}
 	}
 
@@ -282,141 +286,180 @@ double calcParIntQY1X1(double b, double d, double t, double v, double d1,
 
 }
 
-double calcParIntQY1(double b, double d, double t, double v, double d1, double d2, double d3) {
+double calcParIntQY1(double b, double d, double t, double v, double d1,
+		double d2, double d3) {
 	//Fall 4
 	return apply0Int4(F_par, b, d, t, v, d1, d2, d3); ///ACHTUNG noch Falsch
 
 }
 
-double calcOrtIntA(double b, double d, double t, double v, double d1, double d2, double d3) {
+double calcOrtIntA(double b, double d, double t, double v, double d1,
+		double d2, double d3) {
 	return apply0Int4(F_ort, b, d, t, v, d1, d2, d3);
 
 }
 
-double calcOrtIntQY1Y2(double b, double d, double t, double v, double d1, double d2, double d3) {
+double calcOrtIntQY1Y2(double b, double d, double t, double v, double d1,
+		double d2, double d3) {
 	return 0;
 
 }
 
+double calcParIntO0(double b, double d, double t, double v, double d1,
+		double d2, double d3, double eta) {
+	return calcParIntA(b,d,t,v,d1,d2,d3);
+}
+double calcOrtIntO0(double b, double d, double t, double v, double d1,
+		double d2, double d3, double eta) {
+	return calcOrtIntA(b,d,t,v,d1,d2,d3);
+}
 
-double calcParIntO1(double b, double d, double t, double v, double d1, double d2, double d3) {
+double calcParIntO1(double b, double d, double t, double v, double d1,
+		double d2, double d3, double eta) {
 
 	//Elmente anordnen
-	if(b*b+d*d>t*t+v*v){
+	if (b * b + d * d > t * t + v * v) {
 		double tmp = 0;
 
-		tmp = b; b = t; t = tmp;
-		tmp = d; d = v; v = tmp;
+		tmp = b;
+		b = t;
+		t = tmp;
+		tmp = d;
+		d = v;
+		v = tmp;
 
 		d1 = -d1;
 		d2 = -d2;
 		d3 = -d3;
 	}
 
-	return calcParIntA( b, d, t, v, d1, d2, d3);
+	return calcParIntA(b, d, t, v, d1, d2, d3);
 }
 
-double calcOrtIntO1(double b, double d, double t, double v, double d1, double d2, double d3) {
+double calcOrtIntO1(double b, double d, double t, double v, double d1,
+		double d2, double d3, double eta) {
 
 	//Elmente anordnen
-	if(b*b+d*d>t*t+v*v){
+	if (b * b + d * d > t * t + v * v) {
 		double tmp = 0;
 
-		tmp = b; b = t; t = tmp;
-		tmp = d; d = v; v = tmp;
+		tmp = b;
+		b = t;
+		t = tmp;
+		tmp = d;
+		d = v;
+		v = tmp;
 
 		d1 = -d1;
 		d2 = -d2;
 		d3 = -d3;
 	}
 
-	return calcOrtIntA( b, d, t, v, d1, d2, d3);
+	return calcOrtIntA(b, d, t, v, d1, d2, d3);
 }
 
-
-double calcParIntO2(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {
+double calcParIntO2(double b, double d, double t, double v, double d1,
+		double d2, double d3, double eta) {
 
 	//Elmente anordnen
-	if(b*b+d*d>t*t+v*v){
+	if (b * b + d * d > t * t + v * v) {
 		double tmp = 0;
 
-		tmp = b; b = t; t = tmp;
-		tmp = d; d = v; v = tmp;
+		tmp = b;
+		b = t;
+		t = tmp;
+		tmp = d;
+		d = v;
+		v = tmp;
 
 		d1 = -d1;
 		d2 = -d2;
 		d3 = -d3;
 	}
 
-	if((b*b+d*d) *eta <= d1*d1+d2*d2+d3*d3){
-		return calcParIntQX1X2( b, d, t, v, d1, d2, d3);
-	}else{
-		return calcParIntA( b, d, t, v, d1, d2, d3);
+	if ((b * b + d * d) * eta <= d1 * d1 + d2 * d2 + d3 * d3) {
+		return calcParIntQX1X2(b, d, t, v, d1, d2, d3);
+	} else {
+		return calcParIntA(b, d, t, v, d1, d2, d3);
 	}
 
-
 }
 
-double calcOrtIntO2(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {
+double calcOrtIntO2(double b, double d, double t, double v, double d1,
+		double d2, double d3, double eta) {
 
 	//Elmente anordnen
-	if(b*b+d*d>t*t+v*v){
+	if (b * b + d * d > t * t + v * v) {
 		double tmp = 0;
 
-		tmp = b; b = t; t = tmp;
-		tmp = d; d = v; v = tmp;
+		tmp = b;
+		b = t;
+		t = tmp;
+		tmp = d;
+		d = v;
+		v = tmp;
 
 		d1 = -d1;
 		d2 = -d2;
 		d3 = -d3;
 	}
 
-	if((b*b+d*d) *eta <= d1*d1+d2*d2+d3*d3){
+	if ((b * b + d * d) * eta <= d1 * d1 + d2 * d2 + d3 * d3) {
 		return 0;
-	}else{
-		return calcOrtIntA( b, d, t, v, d1, d2, d3);
+	} else {
+		return calcOrtIntA(b, d, t, v, d1, d2, d3);
 	}
 }
 
-
-double calcParIntO3(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {
+double calcParIntO3(double b, double d, double t, double v, double d1,
+		double d2, double d3, double eta) {
 
 	//Achsen vertauschen
-	if(b*b+t*t>d*d+v*v){
+	if (b * b + t * t > d * d + v * v) {
 		double tmp = 0;
 
-		tmp = b; b = d; d = tmp;
-		tmp = t; t = v; v = tmp;
+		tmp = b;
+		b = d;
+		d = tmp;
+		tmp = t;
+		t = v;
+		v = tmp;
 
-		tmp = d1; d1 = d2; d2 = tmp;
+		tmp = d1;
+		d1 = d2;
+		d2 = tmp;
 	}
 
-	if(max(b,t) *eta <= d1){
+	if (max(b, t) * eta <= d1) {
 		return calcParIntQY1X1(b, d, t, v, d1, d2, d3);
-	}else{
-		return calcParIntA( b, d, t, v, d1, d2, d3);
+	} else {
+		return calcParIntA(b, d, t, v, d1, d2, d3);
 	}
 
-
 }
 
-double calcOrtIntO3(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {
+double calcOrtIntO3(double b, double d, double t, double v, double d1,
+		double d2, double d3, double eta) {
 
 	//Elmente anordnen
-	if(b*b+d*d>t*t+v*v){
+	if (b * b + d * d > t * t + v * v) {
 		double tmp = 0;
 
-		tmp = b; b = t; t = tmp;
-		tmp = d; d = v; v = tmp;
+		tmp = b;
+		b = t;
+		t = tmp;
+		tmp = d;
+		d = v;
+		v = tmp;
 
 		d1 = -d1;
 		d2 = -d2;
 		d3 = -d3;
 	}
 
-	if(max(b,t) *eta > d1){
+	if (max(b, t) * eta > d1) {
 		return 0;
-	}else{
-		return calcOrtIntA( b, d, t, v, d1, d2, d3);
+	} else {
+		return calcOrtIntA(b, d, t, v, d1, d2, d3);
 	}
 }
diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp
index f81dfd5..4c63dc4 100644
--- a/src/slpRectangle.hpp
+++ b/src/slpRectangle.hpp
@@ -1,16 +1,11 @@
 #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 g_AY(double, double, double, double);
 // sol = g0(p,y,x,l);
 double g_QY(double, double, double, double);
 
-
 // sol = G00(p,y1,y2,x1,x2,l);
 double G_AY1Y2(double, double, double, double, double, double);
 // sol = G00(p,y1,y2,x1,x2,l);
@@ -18,13 +13,11 @@ double G_AY2X2(double, double, double, double, double, double);
 // sol = G00(p,y1,y2,x1,x2,l);
 double G_QY1Y2(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 = quad0Int4((F_par/F_ort),b,d,t,v,d1,d2,d3);
 double apply0Int4(
 		double(*)(double, double, double, double, double, double, double),
@@ -34,27 +27,37 @@ double apply0Int4(
 double calcParIntA(double, double, double, double, double, double, double);
 double calcParIntQX1X2(double, double, double, double, double, double, double);
 double calcParIntQY1X1(double, double, double, double, double, double, double);
-
 double calcParIntQX1(double, double, double, double, double, double, double);
 
-
 // sol = calcOrtInt.(b,d,t,v,d1,d2,d3);
 double calcOrtIntA(double, double, double, double, double, double, double);
 double calcOrtIntQX1X2(double, double, double, double, double, double, double);
 
+//Voll Analytische Berechnung
+double calcParIntO0(double, double, double, double, double, double, double,
+		double);
+double calcOrtIntO0(double, double, double, double, double, double, double,
+		double);
+
 // Elemente Vertauschen
 // sol = calc...Int01(b,d,t,v,d1,d2,d3);
-double calcParIntO1(double, double, double, double, double, double, double);
-double calcOrtIntO1(double, double, double, double, double, double, double);
+double calcParIntO1(double, double, double, double, double, double, double,
+		double);
+double calcOrtIntO1(double, double, double, double, double, double, double,
+		double);
 
 // Quadratur ueber kleineres Element
 // sol = calc...Int02(b,d,t,v,d1,d2,d3,eta);
-double calcParIntO2(double, double, double, double, double, double, double, double);
-double calcOrtIntO2(double, double, double, double, double, double, double, double);
+double calcParIntO2(double, double, double, double, double, double, double,
+		double);
+double calcOrtIntO2(double, double, double, double, double, double, double,
+		double);
 
 // Quadratur kuerzere Achse/Seite
 // sol = calc...Int03(b,d,t,v,d1,d2,d3,eta);
-double calcParIntO3(double, double, double, double, double, double, double, double);
-double calcOrtIntO3(double, double, double, double, double, double, double, double);
+double calcParIntO3(double, double, double, double, double, double, double,
+		double);
+double calcOrtIntO3(double, double, double, double, double, double, double,
+		double);
 
 #endif
diff --git a/src/test_Fort.m b/src/test_Fort.m
deleted file mode 100644
index ad326d2..0000000
--- a/src/test_Fort.m
+++ /dev/null
@@ -1,19 +0,0 @@
-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_Fpar.m b/src/test_Fpar.m
deleted file mode 100644
index 0f394ef..0000000
--- a/src/test_Fpar.m
+++ /dev/null
@@ -1,44 +0,0 @@
-function [ret] = test_Fpar(a,b,c,d,r,t,u,v,l1,l2,l3)
-%
-% [ret] = test_Fpar(a,b,c,d,r,t,u,v,l1,l2,l3)
-
-    quad_sol = surfDoubleQuad(@(x1,x2,y1,y2)1/sqrt((x1-y1-l1).^2+(x2-y2-l2).^2+l3.^2)...
-                ,a,b,c,d,r,t,u,v,5);
-
-    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);
-      
-      
-      
-      
-      [ mex_Fpar(b,d,t,v,l1,l2,l3) mex_FparY1Y2_2(b,d,t,v,l1,l2,l3);...
-        mex_Fpar(b,d,t,u,l1,l2,l3) mex_FparY1Y2_2(b,d,t,u,l1,l2,l3);...
-        mex_Fpar(b,d,r,v,l1,l2,l3) mex_FparY1Y2_2(b,d,r,v,l1,l2,l3)];
-      
-      [i j] = mex_G00(-0.5, b, d, t + l1, v + l2, l3);
-            
-    sol = intF(@(x1,x2,y1,y2)mex_Fpar(x1,x2,y1,y2,l1,l2,l3),a,b,c,d,r,t,u,v);
-    
-   sol2 = intF(@(x1,x2,y1,y2)mex_FparY1Y2_2(x1,x2,y1,y2,l1,l2,l3),a,b,c,d,r,t,u,v);
-   
-    int2 = @(x1,y1,c,d,u,v)...
-      mex_FparQY2X2(x1,d,y1,v,l1,l2,l3)...
-      -mex_FparQY2X2(x1,c,y1,v,l1,l2,l3)...
-      -mex_FparQY2X2(x1,d,y1,u,l1,l2,l3)...
-      +mex_FparQY2X2(x1,c,y1,u,l1,l2,l3);
-    
-    sol3 = surfQuad(@(x1,y1)(int2(x1,y1,c,d,u,v)),a,b,r,t,5);
-    
-    sol4 =...
-      surfQuad(@(x1,y1)mex_FparQY2X2(x1,d,y1,v,l1,l2,l3),a,b,r,t,5)...
-      -surfQuad(@(x1,y1)mex_FparQY2X2(x1,c,y1,v,l1,l2,l3),a,b,r,t,5)...
-      -surfQuad(@(x1,y1)mex_FparQY2X2(x1,d,y1,u,l1,l2,l3),a,b,r,t,5)...
-      +surfQuad(@(x1,y1)mex_FparQY2X2(x1,c,y1,u,l1,l2,l3),a,b,r,t,5);
-
-    
-    ret = [quad_sol sol sol2 sol3 sol4];
-end
-
diff --git a/src/test_G00.m b/src/test_G00.m
deleted file mode 100644
index 097f8c4..0000000
--- a/src/test_G00.m
+++ /dev/null
@@ -1,19 +0,0 @@
-function [ret] = test_G00(p, x1, x2,l,a,b,c,d)
-%
-% [ret] = test_G00(p, x1, x2,l,a,b,c,d)
-
-    quad_sol = surfQuad(@(y1,y2)((x1-y1).^2+(x2-y2).^2+l.^2).^p,a,b,c,d,4);
-
-
-    [my(1) mq(1)] = mex_G00(p,a,c,x1,x2,l)
-    [my(2) mq(2)] = mex_G00(p,b,d,x1,x2,l)
-    [my(3) mq(3)] = mex_G00(p,b,c,x1,x2,l)
-    [my(4) mq(4)] = mex_G00(p,a,d,x1,x2,l)
-
-
-    my_sol = my(1)+my(2)-my(3)-my(4);
-    mq_sol = mq(1)+mq(2)-mq(3)-mq(4);
-
-    ret = [quad_sol my_sol mq_sol];
-end
-
diff --git a/src/test_SolveSpeed1.fig b/src/test_SolveSpeed1.fig
deleted file mode 100644
index 24401f8..0000000
Binary files a/src/test_SolveSpeed1.fig and /dev/null differ
diff --git a/src/test_functions.m b/src/test_functions.m
deleted file mode 100644
index ce33b3a..0000000
--- a/src/test_functions.m
+++ /dev/null
@@ -1,66 +0,0 @@
-function done=test_functions(eps)
-
-    done = 0;
-%% g0 testen  
-p_g0 = [0.5 0 -0.5 -1 -1.5];
-for l= 1:5
-    for p = p_g0
-        s = test_g0(p,1.4,l,1,4);
-        if(abs(s(1)-s(2))>eps)
-            disp g0
-            l
-            p
-            s
-            return;
-        end
-    end
-end
-    % x muss ausserhalb des Intervalls liegen
-    for p = p_g0
-        s = test_g0(p,1.4,0,2,4);
-        if(abs(s(1)-s(2))>eps)
-            disp g0
-            p
-            s
-            return;
-        end
-    end
-
-%% G00 testen
-p_G00 = [-0.5 -1.5];
-
-    for p = p_G00
-        s = test_G00(p,1.4,0.7,1,2,4,1,3);
-        if(abs(s(1)-s(2))>eps)
-            disp G00
-            p
-            s
-            return;
-        end
-    end
-    
-    for p = p_G00
-        s = test_G00(p,1.1,0.3,0,2,4.5,1,2.5);
-        if(abs(s(1)-s(2))>eps)
-            disp G00
-            p
-            s
-            return;
-        end
-    end
-    
-%% F_par testen
-    
-        
-        s = test_Fpar(1,3,2,3,4,5,1,3,1,3,1);
-        if(abs(s(1)-s(2))>eps)
-            disp Fpar
-            p
-            s
-            return;
-        end
-        
-        
-        
-    done = 1;
-end
\ No newline at end of file
diff --git a/src/test_g0.m b/src/test_g0.m
deleted file mode 100644
index 4163f95..0000000
--- a/src/test_g0.m
+++ /dev/null
@@ -1,26 +0,0 @@
-function [ret] = test_g0(p, x,l,a,b)
-
-
-% quad_sol = quad(@(y)((x-y).^2+l.^2).^p,a,b);
-[no we] = gauss(5,a,b);
-f = @(y)((x-y).^2+l.^2).^p;
-quad_sol = we*f(no);
-
-% cumsum(f(no))'
-
-% X = 0:0.2:5;
-% Y = ((x-X).^2+l.^2).^p;
-% 
-% plot(X,Y)
-
-[my(1) mq(1)] = mex_g0(p,a,x,l);
-[my(2) mq(2)]= mex_g0(p,b,x,l);
-
-my_sol = my(2) - my(1);
-mq_sol = mq(2) - mq(1);
-
-
-
-ret = [quad_sol my_sol mq_sol];
-end
-
diff --git a/src/test_solveErrorS1.m b/src/test_solveErrorS1.m
deleted file mode 100755
index b165322..0000000
--- a/src/test_solveErrorS1.m
+++ /dev/null
@@ -1,114 +0,0 @@
-
-mex mex_build_As1.cpp slpRectangle.cpp
-
-dataIso =[];
-dataAniso =[];
-maxSize = 10^3;
-Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
-% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren
-% Netz = 'exmpl_2DQuad';
-
-%% Isotrope Verfeinerung
-disp 'Isotrop'
-load(Netz)
-
-ref = 0;
-while size(elements,1)<maxSize
-  tic
-  ref = ref+1;
-
-  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
-  A_fine = mex_build_As1(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;
-
-  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
-  marked = 2*ones(1,size(elements,1));
-
-  dataIso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
-  [coordinates, elements] = refineQuad(coordinates,elements,marked);
-  toc
-end
-dataIso
-
-
-%% Anisotrope Verfeinerung
-disp 'Anisotrop'
-load(Netz)
-
-ref = 0;
-while size(elements,1)<maxSize
-  tic
-  ref = ref+1;
-
-  %Netz Isotrop verfeinern (4 Teile)
-  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
-  
-  %Matrix mit feinem Netz berechnen
-  A_fine = mex_build_As1(coordinates_fine,elements_fine);
-%   A2 = mex_build_A(coordinates_fine,elements_fine);
-  
-%   figure(1)
-%   spy(A_fine-A_fine')
-%   
-%   figure (2)
-%   spy(A2-A2')
-%   eh = A_fine - A2;
-%   sum(sum(abs(A2-A2')))
-%   sum(sum(abs(A_fine-A_fine')))
-  
-%   spy(eh)
- 
-%   sum(sum(abs(eh)))
-%   if(eh.^2)
-%     disp 'hä'
-%   end
-  
-  %rechte Seite Aufbauen
-  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
-  
-  %Lösungsvektor bestimmen
-  x_fine = A_fine\b;
-  
-  %Energienorm^2 bestimmen
-  xe_fine = x_fine'*A_fine*x_fine;
-
-  % Welche Elemente sollen verfeinert werden
-  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
-  
-  % Wie soll verfeinert werden
-  marked = mark(x_fine(f2s)',ind,0.5,0.5);
-
-  %Daten zur Auswertung zwischenspeichern (testweise?)
-  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
-  
-  %Sicherheit falls keine Elemente markiert sind
-  if(sum(marked~=1)==0)
-    disp 'Abgebrochen weil nicht mehr verfeinert wird!'
-    break;
-  end
-  
-  %Verfeinern
-  [coordinates, elements] = refineQuad(coordinates,elements,marked);
-  toc
-end
-dataAniso
-
-clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine
-
-
-%% Alles Zeichnen
-figure (1)
-loglog(dataIso(:,1),dataIso(:,2),dataAniso(:,1),dataAniso(:,2));
-
-legend('Isotrop','Anisotrop')
-xlabel('Elemente')
-ylabel('Fehler')
-
-figure (2)
-loglog(dataIso(:,1),dataIso(:,3),dataAniso(:,1),dataAniso(:,3));
-
-legend('Isotrop','Anisotrop')
-xlabel('Elemente')
-ylabel('EnergieNorm')
diff --git a/src/test_solveError_1050_2DQuad.fig b/src/test_solveError_1050_2DQuad.fig
deleted file mode 100644
index 897feaf..0000000
Binary files a/src/test_solveError_1050_2DQuad.fig and /dev/null differ
diff --git a/src/test_solveError_1050_2DQuad.mat b/src/test_solveError_1050_2DQuad.mat
deleted file mode 100644
index abab0a5..0000000
Binary files a/src/test_solveError_1050_2DQuad.mat and /dev/null differ
diff --git a/src/test_solveError_1050_3DFichCube.fig b/src/test_solveError_1050_3DFichCube.fig
deleted file mode 100644
index cf76c24..0000000
Binary files a/src/test_solveError_1050_3DFichCube.fig and /dev/null differ
diff --git a/src/test_solveError_1050_3DFichCube.mat b/src/test_solveError_1050_3DFichCube.mat
deleted file mode 100644
index f4fc430..0000000
Binary files a/src/test_solveError_1050_3DFichCube.mat and /dev/null differ
diff --git a/src/test_solveError_Grid2.m b/src/test_solveError_Grid2.m
deleted file mode 100755
index 52aa416..0000000
--- a/src/test_solveError_Grid2.m
+++ /dev/null
@@ -1,132 +0,0 @@
-
-mex mex_build_As2.cpp slpRectangle.cpp
-
-dataIso =[];
-dataAniso =[];
-maxSize = 500;
-fileName = './test_save';
-Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
-% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren
-% Netz = 'exmpl_2DQuad';
-
-
-%% Anisotrope Verfeinerung
-
-if(~exist('d1'))
-  
-  disp 'Normal'
-  load(Netz)
-
-  ref = 0;
-
-while size(elements,1)<maxSize
-  tic
-  ref = ref+1; 
-  
-
-
-  %Netz Isotrop verfeinern (4 Teile)
-  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
-  
-  %Matrix mit feinem Netz berechnen
-  A_fine = mex_build_A(coordinates_fine,elements_fine);
-  
-  %rechte Seite Aufbauen
-  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
-  
-  %Lösungsvektor bestimmen
-  x_fine = A_fine\b;
-  
-  %Energienorm^2 bestimmen
-  xe_fine = x_fine'*A_fine*x_fine;
-
-  % Welche Elemente sollen verfeinert werden
-  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
-  
-  % Wie soll verfeinert werden
-  marked = mark(x_fine(f2s)',ind,0.5,0.5);
-
-  %Daten zur Auswertung zwischenspeichern (testweise?)
-  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
-  
-  %Sicherheit falls keine Elemente markiert sind
-  if(sum(marked~=1)==0)
-    disp 'Abgebrochen weil nicht mehr verfeinert wird!'
-    break;
-  end
-  
-  %Verfeinern
-  %[coordinates, elements] = refineQuad(coordinates,elements,marked);
-  file = ['save' int2str(ref)];
-  load(file);
-  toc
-end
-d1 = dataAniso
-end
-
-%% Anisotrope Verfeinerung
-disp 'Optimiert'
-load(Netz)
-
-ref = 0;
-while size(elements,1)<maxSize
-  tic
-  ref = ref+1; 
-  
-
-
-  %Netz Isotrop verfeinern (4 Teile)
-  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
-  
-  %Matrix mit feinem Netz berechnen
-  A_fine = mex_build_As2(coordinates_fine,elements_fine);
-  
-  %rechte Seite Aufbauen
-  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
-  
-  %Lösungsvektor bestimmen
-  x_fine = A_fine\b;
-  
-  %Energienorm^2 bestimmen
-  xe_fine = x_fine'*A_fine*x_fine;
-
-  % Welche Elemente sollen verfeinert werden
-  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
-  
-  % Wie soll verfeinert werden
-  marked = mark(x_fine(f2s)',ind,0.5,0.5);
-
-  %Daten zur Auswertung zwischenspeichern (testweise?)
-  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
-  
-  %Sicherheit falls keine Elemente markiert sind
-%   if(sum(marked~=1)==0)
-%     disp 'Abgebrochen weil nicht mehr verfeinert wird!'
-%     break;
-%   end
-  
-  %Verfeinern
-  %[coordinates, elements] = refineQuad(coordinates,elements,marked);
-  file = ['save' int2str(ref)];
-  load(file);
-  toc
-end
-d2 = dataAniso
-
-clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine
-
-
-%% Alles Zeichnen
-figure (1)
-loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));
-
-legend('Normal','optimiert')
-xlabel('Elemente')
-ylabel('Fehler')
-
-figure (2)
-loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));
-
-legend('Normal','Optimiert')
-xlabel('Elemente')
-ylabel('EnergieNorm')
diff --git a/src/test_solveError_Grid3.m b/src/test_solveError_Grid3.m
deleted file mode 100755
index df116a6..0000000
--- a/src/test_solveError_Grid3.m
+++ /dev/null
@@ -1,132 +0,0 @@
-
-mex mex_build_As3.cpp slpRectangle.cpp
-
-dataIso =[];
-dataAniso =[];
-maxSize = 500;
-fileName = './test_save';
-Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
-% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren
-% Netz = 'exmpl_2DQuad';
-
-
-%% Anisotrope Verfeinerung
-
-if(~exist('d1'))
-  
-  disp 'Normal'
-  load(Netz)
-
-  ref = 0;
-
-while size(elements,1)<maxSize
-  tic
-  ref = ref+1; 
-  
-
-
-  %Netz Isotrop verfeinern (4 Teile)
-  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
-  
-  %Matrix mit feinem Netz berechnen
-  A_fine = mex_build_A(coordinates_fine,elements_fine);
-  
-  %rechte Seite Aufbauen
-  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
-  
-  %Lösungsvektor bestimmen
-  x_fine = A_fine\b;
-  
-  %Energienorm^2 bestimmen
-  xe_fine = x_fine'*A_fine*x_fine;
-
-  % Welche Elemente sollen verfeinert werden
-  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
-  
-  % Wie soll verfeinert werden
-  marked = mark(x_fine(f2s)',ind,0.5,0.5);
-
-  %Daten zur Auswertung zwischenspeichern (testweise?)
-  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
-  
-  %Sicherheit falls keine Elemente markiert sind
-  if(sum(marked~=1)==0)
-    disp 'Abgebrochen weil nicht mehr verfeinert wird!'
-    break;
-  end
-  
-  %Verfeinern
-  %[coordinates, elements] = refineQuad(coordinates,elements,marked);
-  file = ['save' int2str(ref)];
-  load(file);
-  toc
-end
-d1 = dataAniso
-end
-
-%% Anisotrope Verfeinerung
-disp 'Optimiert'
-load(Netz)
-
-ref = 0;
-while size(elements,1)<maxSize
-  tic
-  ref = ref+1; 
-  
-
-
-  %Netz Isotrop verfeinern (4 Teile)
-  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
-  
-  %Matrix mit feinem Netz berechnen
-  A_fine = mex_build_As3(coordinates_fine,elements_fine);
-  
-  %rechte Seite Aufbauen
-  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
-  
-  %Lösungsvektor bestimmen
-  x_fine = A_fine\b;
-  
-  %Energienorm^2 bestimmen
-  xe_fine = x_fine'*A_fine*x_fine;
-
-  % Welche Elemente sollen verfeinert werden
-  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
-  
-  % Wie soll verfeinert werden
-  marked = mark(x_fine(f2s)',ind,0.5,0.5);
-
-  %Daten zur Auswertung zwischenspeichern (testweise?)
-  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
-  
-  %Sicherheit falls keine Elemente markiert sind
-%   if(sum(marked~=1)==0)
-%     disp 'Abgebrochen weil nicht mehr verfeinert wird!'
-%     break;
-%   end
-  
-  %Verfeinern
-  %[coordinates, elements] = refineQuad(coordinates,elements,marked);
-  file = ['save' int2str(ref)];
-  load(file);
-  toc
-end
-d2 = dataAniso
-
-clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine
-
-
-%% Alles Zeichnen
-figure (1)
-loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));
-
-legend('Normal','optimiert')
-xlabel('Elemente')
-ylabel('Fehler')
-
-figure (2)
-loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));
-
-legend('Normal','Optimiert')
-xlabel('Elemente')
-ylabel('EnergieNorm')
diff --git a/src/test_solveError_Grids.m b/src/test_solveError_Grids.m
new file mode 100644
index 0000000..5113db9
--- /dev/null
+++ b/src/test_solveError_Grids.m
@@ -0,0 +1,130 @@
+
+mex mex_build_AU.cpp slpRectangle.cpp
+
+dataIso =[];
+dataAniso =[];
+maxSize = 500;
+fileName = './test_save';
+Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
+% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 16.2265 konvergieren
+% Netz = 'exmpl_2DQuad';
+
+
+%% Anisotrope Verfeinerung
+
+if(~exist('d1'))
+  
+  disp 'Normal'
+  load(Netz)
+
+  ref = 0;
+
+while size(elements,1)<maxSize
+  tic
+  ref = ref+1; 
+  
+  %Netz Isotrop verfeinern (4 Teile)
+  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
+  
+  %Matrix mit feinem Netz berechnen
+  A_fine = mex_build_A(coordinates_fine,elements_fine);
+  
+  %rechte Seite Aufbauen
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+  
+  %Lösungsvektor bestimmen
+  x_fine = A_fine\b;
+  
+  %Energienorm^2 bestimmen
+  xe_fine = x_fine'*A_fine*x_fine;
+
+  % Welche Elemente sollen verfeinert werden
+  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
+  
+  % Wie soll verfeinert werden
+  marked = mark(x_fine(f2s)',ind,0.5,0.5);
+
+  %Daten zur Auswertung zwischenspeichern (testweise?)
+  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
+  
+  %Sicherheit falls keine Elemente markiert sind
+  if(sum(marked~=1)==0)
+    disp 'Abgebrochen weil nicht mehr verfeinert wird!'
+    break;
+  end
+  
+  %Verfeinern
+  %[coordinates, elements] = refineQuad(coordinates,elements,marked);
+  file = ['save' int2str(ref)];
+  load(file);
+  toc
+end
+d1 = dataAniso
+end
+
+%% Anisotrope Verfeinerung
+disp 'Optimiert'
+load(Netz)
+
+ref = 0;
+while size(elements,1)<maxSize
+  tic
+  ref = ref+1; 
+  
+
+
+  %Netz Isotrop verfeinern (4 Teile)
+  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
+  
+  %Matrix mit feinem Netz berechnen
+  A_fine = mex_build_AU(coordinates_fine,elements_fine);
+  
+  %rechte Seite Aufbauen
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+  
+  %Lösungsvektor bestimmen
+  x_fine = A_fine\b;
+  
+  %Energienorm^2 bestimmen
+  xe_fine = x_fine'*A_fine*x_fine;
+
+  % Welche Elemente sollen verfeinert werden
+  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
+  
+  % Wie soll verfeinert werden
+  marked = mark(x_fine(f2s)',ind,0.5,0.5);
+
+  %Daten zur Auswertung zwischenspeichern (testweise?)
+  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
+  
+  %Sicherheit falls keine Elemente markiert sind
+%   if(sum(marked~=1)==0)
+%     disp 'Abgebrochen weil nicht mehr verfeinert wird!'
+%     break;
+%   end
+  
+  %Verfeinern
+  %[coordinates, elements] = refineQuad(coordinates,elements,marked);
+  file = ['save' int2str(ref)];
+  load(file);
+  toc
+end
+d2 = dataAniso
+
+clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine
+
+
+%% Alles Zeichnen
+figure (1)
+loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));
+
+legend('Normal','optimiert')
+xlabel('Elemente')
+ylabel('Fehler')
+
+figure (2)
+loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));
+
+legend('Normal','Optimiert')
+xlabel('Elemente')
+ylabel('EnergieNorm')
diff --git a/src/test_solveS.m b/src/test_solveS.m
deleted file mode 100644
index 96d646e..0000000
--- a/src/test_solveS.m
+++ /dev/null
@@ -1,33 +0,0 @@
-
-N = size(elements,1);
-
-A1 = zeros(N);
-
-x = zeros(N,2);
-xe = zeros(1,2);
-t = zeros(1,2);
-
-% untere schranke s t obere schranke k l
-intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ...
-    f(k1,k2,l1,l2)-f(k1,k2,l1,t2)-f(k1,k2,t1,l2)+f(k1,k2,t1,t2)...
-    -f(k1,s2,l1,l2)+f(k1,s2,l1,t2)+f(k1,s2,t1,l2)-f(k1,s2,t1,t2)...
-    -f(s1,k2,l1,l2)+f(s1,k2,l1,t2)+f(s1,k2,t1,l2)-f(s1,l2,t1,t2)...
-    +f(s1,s2,l1,l2)-f(s1,0,l1,t2)-f(s1,s2,t1,l2)+f(s1,s2,t1,t2);
-tic
-A2 = build_A2(coordinates,elements);
-t(1) = toc;
-disp ' '
-tic
-A2 = mex_build_A(coordinates,elements);
-t(2) = toc;
-b = sqrt(sum(quadNorm(coordinates,elements,'w').^2,2));
-
-x(:,1) = A1\b;
-xe(1) = x(:,1)'*A1*x(:,1);
-
-x(:,2) = A2\b;
-xe(2) = x(:,2)'*A2*x(:,2);
-
-x;
-xe
-t
\ No newline at end of file
diff --git a/src/test_solveSpeed.m b/src/test_solveSpeed.m
deleted file mode 100644
index eadacaa..0000000
--- a/src/test_solveSpeed.m
+++ /dev/null
@@ -1,23 +0,0 @@
-
-rt_times = 5;
-
-time = zeros(rt_times,3);
-
-    test_solveS
-    
-    time(1,:) = [size(x,1) t];
-
-for l = 2:rt_times
-    test_refine
-    test_solveS
-    
-    time(l,:) = [size(x,1) t];
-end
-
-
-loglog(time(:,1),time(:,2),time(:,1),time(:,3));
-
-legend('Matlab','Mex')
-xlabel('Elemente')
-ylabel('Zeit')
-
diff --git a/src/test_solveSpeed1.mat b/src/test_solveSpeed1.mat
deleted file mode 100644
index 904568a..0000000
Binary files a/src/test_solveSpeed1.mat and /dev/null differ