From d6aa09af7498e604491e34617b2cbdcc57d750bd Mon Sep 17 00:00:00 2001
From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Fri, 26 Aug 2011 07:46:48 +0000
Subject: [PATCH] [src] Elementdrehen angefangen -> leider Fehlerhaft

git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@36 26120e32-c555-405d-b3e1-1f783fb42516
---
 src/mex_build_As1.cpp   | 329 ++++++++++++++++++++++++++++++++++++++++
 src/test_solveError.m   |  17 ++-
 src/test_solveErrorS1.m | 103 +++++++++++++
 3 files changed, 448 insertions(+), 1 deletion(-)
 create mode 100644 src/mex_build_As1.cpp
 create mode 100755 src/test_solveErrorS1.m

diff --git a/src/mex_build_As1.cpp b/src/mex_build_As1.cpp
new file mode 100644
index 0000000..ec82bee
--- /dev/null
+++ b/src/mex_build_As1.cpp
@@ -0,0 +1,329 @@
+#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;
+
+	int itmp;
+
+	double * vtmp;
+
+	//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];
+
+
+			if(ry = rx){
+
+				if(xa[rxa]*xa[rxa]+xb[rxb]*xb[rxb]>ya[rya]*ya[rya]+yb[ryb]*yb[ryb]){
+
+	/*				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,d[0],d[1],d[2]);
+					printf(" | ");
+*/
+
+					printf("(%3.d %3.d)",j,k);
+
+					vtmp = xa; xa = ya; ya = vtmp;
+					vtmp = xb; xb = yb; yb = vtmp;
+
+					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];
+
+					/*
+
+					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,d[0],d[1],d[2]);
+					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("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);
+				//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
+				if (rxa == rya) {
+					//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],
+					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
+					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
+					tmp
+							= apply0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(ya[rxa]), fabs(yb[rxb]), d[rxa],
+									d[rxb], d[rx]);
+					//	 printf("%d%d|%.2f\n",j,k,tmp);
+				} else {
+					//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],
+					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
+					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
+					tmp
+							= apply0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(yb[rxa]), fabs(ya[rxb]), d[rxa],
+									d[rxb], d[rx]);
+					//	 printf("%d%d|%.2f\n",j,k,tmp);
+				}
+
+			} else {
+				//printf("(%d,%d@%d,%d)", rxa, rxb, rya, ryb);
+				//printf("(%d,%d)", rx, ry);
+				//printf("\n");
+
+				if (rxa == rya) {
+					tmp
+							= apply0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),
+									fabs(ya[rya]), fabs(yb[ryb]), d[rxb],
+									d[rxa], d[rx]);
+				} else if (rxa == ryb) {
+					tmp
+							= apply0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),
+									fabs(yb[ryb]), fabs(ya[rya]), d[rxb],
+									d[rxa], d[rx]);
+				} else if (rxb == rya) {
+					tmp
+							= apply0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(ya[rya]), fabs(yb[ryb]), d[rxa],
+									d[rxb], d[rx]);
+				} else {
+					tmp
+							= apply0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),
+									fabs(yb[ryb]), fabs(ya[rya]), d[rxa],
+									d[rxb], d[rx]);
+				}
+				//A[(k * em) + j] = tmp;
+
+			}
+			A[(k * em) + j] = 1. / (4 * M_PI) * tmp;
+
+		}
+
+	}
+	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/test_solveError.m b/src/test_solveError.m
index 0316541..fa537bf 100644
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -1,7 +1,7 @@
 
 dataIso =[];
 dataAniso =[];
-maxSize = 1050;
+maxSize = 10^2;
 Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
 % Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren
 % Netz = 'exmpl_2DQuad';
@@ -40,22 +40,37 @@ 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);
   toc
 end
diff --git a/src/test_solveErrorS1.m b/src/test_solveErrorS1.m
new file mode 100755
index 0000000..4143f9c
--- /dev/null
+++ b/src/test_solveErrorS1.m
@@ -0,0 +1,103 @@
+
+dataIso =[];
+dataAniso =[];
+maxSize = 10^2;
+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);
+  eh = A_fine-mex_build_A(coordinates_fine,elements_fine);
+  figure (3)
+  spy(eh)
+ 
+  sum(sum(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')
-- 
2.47.3