From c52899ecb5e0b1a9105de6f8f948ccdc9ae0c0da Mon Sep 17 00:00:00 2001
From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Mon, 29 Aug 2011 10:12:23 +0000
Subject: [PATCH] [src] Elemente Vertauschen geht aber ohne sichtbaren erfolg
 [src] inneres Element durch Quadratur... geht noch nicht

git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@39 26120e32-c555-405d-b3e1-1f783fb42516
---
 src/mex_Fpar.cpp                              |   2 +-
 src/mex_G00.cpp                               |   2 +-
 src/mex_build_A.cpp                           |  20 +-
 src/mex_build_As1.cpp                         | 101 +++---
 src/mex_build_As2.cpp                         | 328 ++++++++++++++++++
 src/mex_g0.cpp                                |   2 +-
 src/mexit.m                                   |   2 +
 src/slpRectangle.cpp                          |  70 +++-
 src/slpRectangle.hpp                          |  20 +-
 src/test_Fpar.m                               |  22 +-
 src/test_G00.m                                |  16 +-
 src/test_g0.m                                 |  12 +-
 src/test_solveErrorS1.m                       |  16 +-
 ...veError_Grid.m => test_solveError_CGrid.m} |   0
 src/test_solveError_Grid1.m                   | 130 +++++++
 src/test_solveError_Grid2.m                   | 116 ++++---
 16 files changed, 719 insertions(+), 140 deletions(-)
 create mode 100644 src/mex_build_As2.cpp
 rename src/{test_solveError_Grid.m => test_solveError_CGrid.m} (100%)
 create mode 100755 src/test_solveError_Grid1.m

diff --git a/src/mex_Fpar.cpp b/src/mex_Fpar.cpp
index c02cfb6..9103e4d 100644
--- a/src/mex_Fpar.cpp
+++ b/src/mex_Fpar.cpp
@@ -31,7 +31,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	//sol = G00(p,y1,y2,x1,x2,l);
 
 	*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]));
-	*mxGetPr(plhs[1]) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+	*mxGetPr(plhs[1]) = F_parY1Y2_2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
 
 //	printf("%f",FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6])));
 
diff --git a/src/mex_G00.cpp b/src/mex_G00.cpp
index 9f9fd98..24379f5 100644
--- a/src/mex_G00.cpp
+++ b/src/mex_G00.cpp
@@ -29,7 +29,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	//sol = G00(p,y1,y2,x1,x2,l);
 
 	*mxGetPr(plhs[0]) = G00(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
-	*mxGetPr(plhs[1]) = slpADLO(*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
+	*mxGetPr(plhs[1]) = doubleQuad(*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
index 48a30b0..0cdec5f 100644
--- a/src/mex_build_A.cpp
+++ b/src/mex_build_A.cpp
@@ -1,3 +1,11 @@
+/*
+ * voll Analytisch
+ * ohne Optimierungen
+ *
+ *
+ */
+
+
 #include <iostream>
 #include <cmath>
 #include <cassert>
@@ -225,7 +233,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					//		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]),
+							= apply0Int4(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);
@@ -234,7 +242,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					//		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]),
+							= apply0Int4(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);
@@ -247,22 +255,22 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 				if (rxa == rya) {
 					tmp
-							= apply0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),
+							= apply0Int4(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]),
+							= apply0Int4(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]),
+							= apply0Int4(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]),
+							= apply0Int4(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),
 									fabs(yb[ryb]), fabs(ya[rya]), d[rxa],
 									d[rxb], d[rx]);
 				}
diff --git a/src/mex_build_As1.cpp b/src/mex_build_As1.cpp
index 485b74f..2356c72 100644
--- a/src/mex_build_As1.cpp
+++ b/src/mex_build_As1.cpp
@@ -1,3 +1,11 @@
+/*
+ * voll Analytisch
+ * Element mit groesserer Diagonale nach aussen gedreht
+ *
+ *
+ */
+
+
 #include <iostream>
 #include <cmath>
 #include <cassert>
@@ -72,11 +80,14 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	double tmp,tmp2;
 
 	int itmp;
-
+	double dtmp;
 	double * vtmp;
 
+	long double xda,xdb,yda,ydb;
+
 	//LageInformationen
 	int rx, rxa, rxb, ry, rya, ryb;
+	int rtxa, rtxb, rtx;
 
 	//Ausrechnen
 	for (j = 0; j < em; ++j) {
@@ -110,17 +121,17 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 		// Lageeigenschaften des Flaechenstuecks
 		//if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])
 
-		rx = dimOfVec(xn);
+		rtx = dimOfVec(xn);
 
-		rxa = dimOfVec(xa);
+		rtxa = dimOfVec(xa);
 
-		rxb = dimOfVec(xb);
+		rtxb = dimOfVec(xb);
 
 
 		//kleinste Ecke finden und fuer \delta verwenden
 
-		if (xa[rxa] > 0) {
-			if (xb[rxb] > 0) {
+		if (xa[rtxa] > 0) {
+			if (xb[rtxb] > 0) {
 				for (i = 0; i < 3; ++i)
 					dt[i] = -x0[i];
 			} else {
@@ -128,7 +139,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					dt[i] = -x3[i];
 			}
 		} else {
-			if (xb[rxb] > 0) {
+			if (xb[rtxb] > 0) {
 				for (i = 0; i < 3; ++i)
 					dt[i] = -x1[i];
 			} else {
@@ -145,6 +156,15 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 		// printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);
 
 		for (k = 0; k < em; ++k) {
+
+			xda = xa[rtxa];
+			xdb = xb[rtxb];
+
+			rx = rtx;
+			rxa = rtxa;
+			rxb = rtxb;
+
+
 			y0[0] = C[(int) E[k] - 1];
 			y0[1] = C[cm + (int) E[k] - 1];
 			y0[2] = C[2 * cm + (int) E[k] - 1];
@@ -198,22 +218,20 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 				}
 			}
 
-			//	for (i = 0; i<3;++i)
-			//		d[i] = y0[i]-x0[i];
-
+			yda = ya[rya];
+			ydb = yb[ryb];
 
-			if(ry == rx){
 
-				if(xa[rxa]*xa[rxa]+xb[rxb]*xb[rxb]>ya[rya]*ya[rya]+yb[ryb]*yb[ryb]){
+				if(xda*xda+xdb*xdb>yda*yda+ydb*ydb){
 
-					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@%d,%d)",rxa,rxb,rya,ryb);
+					printf("[%.2f,%.2f@%.2f,%.2f]",xda,xb[rxb],ya[rya],yb[ryb]);
 					printf("{%d%d} %.1f %.1f %.1f",j,k,d[0],d[1],d[2]);
-					printf("\n");
+					printf("\n");*/
 
 
-					vtmp = xa; xa = ya; ya = vtmp;
-					vtmp = xb; xb = yb; yb = vtmp;
+					dtmp = xda; xda = yda; yda = dtmp;
+					dtmp = xdb; xdb = ydb; ydb = dtmp;
 
 					itmp = rxa; rxa = rya; rya = itmp;
 					itmp = rxb; rxb = ryb; ryb = itmp;
@@ -224,19 +242,17 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 
 
-					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@%d,%d)",rxa,rxb,rya,ryb);
+					printf("[%.2f,%.2f@%.2f,%.2f]",xda,xb[rxb],ya[rya],yb[ryb]);
 					printf("{%d%d} %.1f %.1f %.1f",j,k,d[0],d[1],d[2]);
-					printf("\n");
+					printf("\n");*/
 				}
 
 
-			}
-
 			//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("[%.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) {
@@ -244,20 +260,20 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					//		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]);
 
-					if(xa[rxa]*xa[rxa]+xb[rxb]*xb[rxb]<=ya[rya]*ya[rya]+yb[ryb]*yb[ryb]){
+			//		if(xda*xda+xdb*xdb<=yda*yda+ydb*ydb){
 
 					tmp
-							= apply0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),
-									fabs(ya[rya]), fabs(yb[ryb]), d[rxa],
+							= apply0Int4(F_par, fabs(xda), fabs(xdb),
+									fabs(yda), fabs(ydb), d[rxa],
 									d[rxb], d[rx]);
 
-					}else{
+			//		}else{
 
-					tmp
-							= apply0Int(F_par, fabs(ya[rya]), fabs(yb[ryb]),
-									fabs(xa[rxa]), fabs(xb[rxb]), -d[rxa],
-									-d[rxb], -d[rx]);
-					}
+			//		tmp
+			//				= apply0Int4(F_par, fabs(yda), fabs(ydb),
+			//						fabs(xda), fabs(xdb), -d[rxa],
+			//						-d[rxb], -d[rx]);
+			//		}
 			//		if(fabs(tmp-tmp2)>10e-16)
 			//			printf("wtf");
 					//	 printf("%d%d|%.2f\n",j,k,tmp);
@@ -266,8 +282,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					//		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],
+							= apply0Int4(F_par, fabs(xda), fabs(xdb),
+									fabs(ydb), fabs(yda), d[rxa],
 									d[rxb], d[rx]);
 					//	 printf("%d%d|%.2f\n",j,k,tmp);
 				}
@@ -279,26 +295,25 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 				if (rxa == rya) {
 					tmp
-							= apply0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),
-									fabs(ya[rya]), fabs(yb[ryb]), d[rxb],
+							= apply0Int4(F_ort, fabs(xdb), fabs(xda),
+									fabs(yda), fabs(ydb), 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],
+							= apply0Int4(F_ort, fabs(xdb), fabs(xda),
+									fabs(ydb), fabs(yda), 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],
+							= apply0Int4(F_ort, fabs(xda), fabs(xdb),
+									fabs(yda), fabs(ydb), 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],
+							= apply0Int4(F_ort, fabs(xda), fabs(xdb),
+									fabs(ydb), fabs(yda), d[rxa],
 									d[rxb], d[rx]);
 				}
-				//A[(k * em) + j] = tmp;
 
 			}
 			A[(k * em) + j] = 1. / (4 * M_PI) * tmp;
diff --git a/src/mex_build_As2.cpp b/src/mex_build_As2.cpp
new file mode 100644
index 0000000..c602509
--- /dev/null
+++ b/src/mex_build_As2.cpp
@@ -0,0 +1,328 @@
+/*
+ * 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;
+	int rtxa, rtxb, rtx;
+
+	//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])
+
+		rtx = dimOfVec(xn);
+
+		rtxa = dimOfVec(xa);
+
+		rtxb = dimOfVec(xb);
+
+
+		//kleinste Ecke finden und fuer \delta verwenden
+
+		if (xa[rtxa] > 0) {
+			if (xb[rtxb] > 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[rtxb] > 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) {
+
+			xda = xa[rtxa];
+			xdb = xb[rtxb];
+
+			rx = rtx;
+			rxa = rtxa;
+			rxb = rtxb;
+
+
+			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];
+
+				//groesseres Element nach aussen drehen
+				if(xda*xda+xdb*xdb>yda*yda+ydb*ydb){
+
+					dtmp = xda; xda = yda; yda = dtmp;
+					dtmp = xdb; xdb = ydb; ydb = dtmp;
+
+					itmp = rxa; rxa = rya; rya = itmp;
+					itmp = rxb; rxb = ryb; ryb = itmp;
+					itmp = rx; rx = ry; ry = itmp;
+
+					for (i = 0; i<3;++i)
+						d[i] = -d[i];
+
+				}
+
+				//if(xda*xda+xdb*xdb>=d[0]*d[0]+d[1]*d[1]+d[2]*d[2]){
+				//	printf(".");
+				//}
+
+			//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]);
+
+			//		if(xda*xda+xdb*xdb<d[0]*d[0]+d[1]*d[1]+d[2]*d[2]){
+
+					tmp
+							= apply0Int4(F_par, fabs(xda), fabs(xdb),
+									fabs(yda), fabs(ydb), d[rxa],
+									d[rxb], d[rx]);
+
+			//		}else{
+
+					tmp2
+							= apply0Int2(F_parY1Y2_2, fabs(xda), fabs(xdb),
+									fabs(yda), fabs(ydb), d[rxa],
+									d[rxb], d[rx]);
+			//		}
+					if(fabs(tmp-tmp2)>10e-8)
+						return;
+					//	 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
+							= apply0Int4(F_par, 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
+							= apply0Int4(F_ort, fabs(xdb), fabs(xda),
+									fabs(yda), fabs(ydb), d[rxb],
+									d[rxa], d[rx]);
+				} else if (rxa == ryb) {
+					tmp
+							= apply0Int4(F_ort, fabs(xdb), fabs(xda),
+									fabs(ydb), fabs(yda), d[rxb],
+									d[rxa], d[rx]);
+				} else if (rxb == rya) {
+					tmp
+							= apply0Int4(F_ort, fabs(xda), fabs(xdb),
+									fabs(yda), fabs(ydb), d[rxa],
+									d[rxb], d[rx]);
+				} else {
+					tmp
+							= apply0Int4(F_ort, 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_g0.cpp b/src/mex_g0.cpp
index 9924813..7dfb788 100644
--- a/src/mex_g0.cpp
+++ b/src/mex_g0.cpp
@@ -27,7 +27,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	//printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
 
 	*mxGetPr(plhs[0]) = g0(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
-	*mxGetPr(plhs[1]) = compute_g0(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+	*mxGetPr(plhs[1]) = *mxGetPr(prhs[1])!=0?singleQuad(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3])):0;
 
 
 }
diff --git a/src/mexit.m b/src/mexit.m
index ed3eaad..f868508 100644
--- a/src/mexit.m
+++ b/src/mexit.m
@@ -5,4 +5,6 @@ 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/slpRectangle.cpp b/src/slpRectangle.cpp
index 42f480a..b92c6a4 100644
--- a/src/slpRectangle.cpp
+++ b/src/slpRectangle.cpp
@@ -3,12 +3,17 @@
 #include <cassert>
 #include <stdlib.h>
 
+#include <mex.h>
+
 #include "slpRectangle.hpp"
 
-#define EPS 0.00001
+#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};
+
 int sign(double);
 //double arsinh(double);
 
@@ -20,6 +25,31 @@ int inline sign(double x) {
  return log(x + sqrt(x * x + 1));
  }
  */
+
+double singleQuad(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;
+	}
+
+	return sol;
+}
+
+
+double doubleQuad(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];
+		}
+	}
+	return sol;
+}
+
 //y-x muss != 0 sein
 double g0(double p, double y, double x, double l) {
 	//printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);
@@ -155,6 +185,29 @@ double F_parX1Y2(double x1, double x2, double y1, double y2, double d1, double d
 	return sol;
 }
 
+double F_parY1Y2_2(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 *= doubleQuad(-0.5, x1, x2, y1 + d1, y2 + d2, d3);
+
+	if ((x1 - y1 - d1) != 0)
+		sol -= (x1 - y1 - d1) * singleQuad(0.5, x1, y1 + d1,
+				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));
+
+	if ((x2 - y2 - d2) != 0)
+		sol -= (x2 - y2 - d2) * singleQuad(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) {
 
@@ -180,7 +233,7 @@ double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,
 	return sol / 2.;
 }
 
-double applyInt(
+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) {
@@ -205,7 +258,7 @@ double applyInt(
 	 */
 }
 
-double apply0Int(
+double apply0Int4(
 		double(*f)(double, double, double, double, double, double, double),
 		double b, double d, double t, double v, double d1, double d2, double d3) {
 
@@ -219,6 +272,17 @@ double apply0Int(
 
 }
 
+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);
+
+}
+
 double slpADLO(double y1, double y2, double x1, double x2, double a) {
 	double G3 = 0;
 	double gL = 0;
diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp
index 87e3bdf..5b3a715 100644
--- a/src/slpRectangle.hpp
+++ b/src/slpRectangle.hpp
@@ -4,6 +4,12 @@
 int sign(double);
 double arsinh(double);
 
+// sol = g0(p,y,x,l); Quadratur !!! 0-y
+double singleQuad(double, double, double, double);
+
+// sol = G00(p,y1,y2,x1,x2,l); Quadratur !!! 0-y1 & 0-y2
+double doubleQuad(double, double, double, double, double, double);
+
 // sol = g0(p,y,x,l);
 double g0(double, double, double, double);
 // sol = G00(p,y1,y2,x1,x2,l);
@@ -23,19 +29,29 @@ double F_parY1X2(double, double, double, double, double, double, double);
 // sol = F_par(x1,x2,y1,y2,d1,d2,d3);
 double F_parX1Y2(double, double, double, double, double, double, double);
 
+// sol = F_par(x1,x2,y1,y2,d1,d2,d3); Quadratur!!!
+double F_parY1Y2_2(double, double, double, double, double, double, double);
+
+/*
+// funktionen von Maik
 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 applyInt4(
 		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 apply0Int4(
+		double(*)(double, double, double, double, double, double, double),
+		double, double, double, double, double, double, double);
+
+double apply0Int2(
 		double(*)(double, double, double, double, double, double, double),
 		double, double, double, double, double, double, double);
 
diff --git a/src/test_Fpar.m b/src/test_Fpar.m
index 27c3e6c..e1f70a7 100644
--- a/src/test_Fpar.m
+++ b/src/test_Fpar.m
@@ -2,18 +2,32 @@ function [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,4);
+                ,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);
+      
+      
+    intF2 = @(f,a,b,c,d,r,t,u,v)...
+        f(b,d,t,v)-f(a,d,t,v)-f(b,c,t,v)+f(a,c,t,v);
+      
+      
+     [dm sol] = mex_Fpar(b,d,t,v,l1,l2,l3)
+
+     [dm sol1] =  mex_Fpar(b,d,0,v,l1,l2,l3)
+     sol = sol - sol1;
+     [dm sol1] =  mex_Fpar(b,d,t,0,l1,l2,l3)
+     sol = sol - sol1;
+     [dm sol1] =  mex_Fpar(b,d,0,0,l1,l2,l3)
+    q_sol = sol + sol1;
+      
+      
             
     sol = intF(@(x1,x2,y1,y2)mex_Fpar(x1,x2,y1,y2,l1,l2,l3),a,b,c,d,r,t,u,v);
 
-
-
-    ret = [quad_sol sol];
+    ret = [quad_sol sol q_sol sol/q_sol];
 end
 
diff --git a/src/test_G00.m b/src/test_G00.m
index c23354f..66e7f7a 100644
--- a/src/test_G00.m
+++ b/src/test_G00.m
@@ -4,20 +4,14 @@ function [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) mk(1)] = mex_G00(p,a,c,x1,x2,l);
-    [my(2) mk(2)] = mex_G00(p,b,d,x1,x2,l);
-    [my(3) mk(3)] = mex_G00(p,b,c,x1,x2,l);
-    [my(4) mk(4)] = mex_G00(p,a,d,x1,x2,l);
+    [my(1) ] = mex_G00(p,a,c,x1,x2,l);
+    [my(2) mq_sol] = mex_G00(p,b,d,x1,x2,l);
+    [my(3) ] = mex_G00(p,b,c,x1,x2,l);
+    [my(4) ] = mex_G00(p,a,d,x1,x2,l);
 
 
     my_sol = my(1)+my(2)-my(3)-my(4);
-    if(p==-0.5)
-        mk_sol = mk(1)+mk(2)-mk(3)-mk(4);
-    else
-        mk_sol = nan;
-    end
 
-
-    ret = [quad_sol my_sol mk_sol];
+    ret = [quad_sol my_sol mq_sol];
 end
 
diff --git a/src/test_g0.m b/src/test_g0.m
index 89e2a02..dbd75dc 100644
--- a/src/test_g0.m
+++ b/src/test_g0.m
@@ -2,22 +2,24 @@ function [ret] = test_g0(p, x,l,a,b)
 
 
 % quad_sol = quad(@(y)((x-y).^2+l.^2).^p,a,b);
-[no we] = gauss(11,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) mk(1)] = mex_g0(p,a,x,l);
-[my(2) mk(2)] = mex_g0(p,b,x,l);
+my(1) = mex_g0(p,a,x,l);
+[my(2) mq_sol]= mex_g0(p,b,x,l);
 
 my_sol = my(2) - my(1);
-mk_sol = mk(2) - mk(1);
 
 
-ret = [quad_sol my_sol mk_sol];
+
+ret = [quad_sol my_sol mq_sol];
 end
 
diff --git a/src/test_solveErrorS1.m b/src/test_solveErrorS1.m
index c367bb2..b165322 100755
--- a/src/test_solveErrorS1.m
+++ b/src/test_solveErrorS1.m
@@ -3,7 +3,7 @@ mex mex_build_As1.cpp slpRectangle.cpp
 
 dataIso =[];
 dataAniso =[];
-maxSize = 10^2;
+maxSize = 10^3;
 Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
 % Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren
 % Netz = 'exmpl_2DQuad';
@@ -47,23 +47,23 @@ while size(elements,1)<maxSize
   
   %Matrix mit feinem Netz berechnen
   A_fine = mex_build_As1(coordinates_fine,elements_fine);
-  A2 = mex_build_A(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;
+%   eh = A_fine - A2;
 %   sum(sum(abs(A2-A2')))
 %   sum(sum(abs(A_fine-A_fine')))
   
-  spy(eh)
+%   spy(eh)
  
-  sum(sum(abs(eh)))
-  if(eh.^2)
-    disp 'hä'
-  end
+%   sum(sum(abs(eh)))
+%   if(eh.^2)
+%     disp 'hä'
+%   end
   
   %rechte Seite Aufbauen
   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
diff --git a/src/test_solveError_Grid.m b/src/test_solveError_CGrid.m
similarity index 100%
rename from src/test_solveError_Grid.m
rename to src/test_solveError_CGrid.m
diff --git a/src/test_solveError_Grid1.m b/src/test_solveError_Grid1.m
new file mode 100755
index 0000000..fe04315
--- /dev/null
+++ b/src/test_solveError_Grid1.m
@@ -0,0 +1,130 @@
+
+mex mex_build_As1.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_As1(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_Grid2.m b/src/test_solveError_Grid2.m
index f90d5af..52aa416 100755
--- a/src/test_solveError_Grid2.m
+++ b/src/test_solveError_Grid2.m
@@ -1,4 +1,6 @@
 
+mex mex_build_As2.cpp slpRectangle.cpp
+
 dataIso =[];
 dataAniso =[];
 maxSize = 500;
@@ -8,61 +10,15 @@ Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
 % Netz = 'exmpl_2DQuad';
 
 
-% %% Anisotrope Verfeinerung
-% disp 'Anisotrop1'
-% 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
+%% Anisotrope Verfeinerung
 
+if(~exist('d1'))
+  
+  disp 'Normal'
+  load(Netz)
 
-%% Anisotrope Verfeinerung
-disp 'Anisotrop2'
-load(Netz)
+  ref = 0;
 
-ref = 0;
 while size(elements,1)<maxSize
   tic
   ref = ref+1; 
@@ -73,7 +29,7 @@ while size(elements,1)<maxSize
   [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);
+  A_fine = mex_build_A(coordinates_fine,elements_fine);
   
   %rechte Seite Aufbauen
   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
@@ -105,6 +61,56 @@ while size(elements,1)<maxSize
   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
@@ -114,13 +120,13 @@ clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine
 figure (1)
 loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));
 
-legend('Isotrop','Anisotrop')
+legend('Normal','optimiert')
 xlabel('Elemente')
 ylabel('Fehler')
 
 figure (2)
 loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));
 
-legend('Isotrop','Anisotrop')
+legend('Normal','Optimiert')
 xlabel('Elemente')
 ylabel('EnergieNorm')
-- 
2.47.3