From 742a7f8b008beda36c39b62c27afc47f0c49fced Mon Sep 17 00:00:00 2001
From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Thu, 25 Aug 2011 13:47:08 +0000
Subject: [PATCH] =?utf8?q?[src]=20Analytisch=20Y1X2=20eingebaut=20sowie=20?=
 =?utf8?q?X1Y2=20=20->=20f=C3=BCr=20f=C3=BCr=20paralelle=20Elemente?=
MIME-Version: 1.0
Content-Type: text/plain; charset=utf8
Content-Transfer-Encoding: 8bit

git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@35 26120e32-c555-405d-b3e1-1f783fb42516
---
 src/mex_FparY1X2.cpp  | 44 +++++++++++++++++++++++++++++++++++++++
 src/slpRectangle.cpp  | 48 ++++++++++++++++++++++++++++++++++++++++++-
 src/slpRectangle.hpp  | 11 +++++++++-
 src/test_solveError.m | 20 ++++++++++++------
 4 files changed, 115 insertions(+), 8 deletions(-)
 create mode 100755 src/mex_FparY1X2.cpp

diff --git a/src/mex_FparY1X2.cpp b/src/mex_FparY1X2.cpp
new file mode 100755
index 0000000..81453f1
--- /dev/null
+++ b/src/mex_FparY1X2.cpp
@@ -0,0 +1,44 @@
+#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);
+	plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
+
+//	plhs[0] = mxCreateDoubleMatrix(1, 2, mxREAL);
+
+
+	//printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+
+	//sol = G00(p,y1,y2,x1,x2,l);
+
+	*mxGetPr(plhs[0]) = F_parX1Y2(*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]));
+
+//	printf("%f",FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6])));
+
+//	*(mxGetPr(plhs[0])+1) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+
+
+}
+
+
+
diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp
index 7e9a7b0..42f480a 100644
--- a/src/slpRectangle.cpp
+++ b/src/slpRectangle.cpp
@@ -86,7 +86,7 @@ double G00(double p, double y1, double y2, double x1, double x2, double l) {
 	return sol;
 }
 
-double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
+double F_parY1Y2(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);
@@ -109,6 +109,52 @@ double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
 	return sol;
 }
 
+double F_parY1X2(double x1, double x2, double y1, double y2, double d1, double d2,
+		double d3) {
+
+	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);
+	double sol = (x1 - y1 - d1) * (x2 - y2 - d2);
+
+	if (sol != 0)
+		sol *= G00(-0.5, x1, -y2 - d2, y1 + d1, -x2, d3);
+
+	if ((x1 - y1 - d1) != 0)
+		sol -= (x1 - y1 - d1) * g0(0.5, x1, y1 + d1,
+				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));
+
+	if ((x2 - y2 - d2) != 0)
+		sol -= (x2 - y2 - d2) * g0(0.5, -y2 - d2, -x2,
+				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_parX1Y2(double x1, double x2, double y1, double y2, double d1, double d2,
+		double d3) {
+
+	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);
+	double sol = (x1 - y1 - d1) * (x2 - y2 - d2);
+
+	if (sol != 0)
+		sol *= G00(-0.5, -y1 - d1, x2, -x1, y2 + d2, d3);
+
+	if ((x1 - y1 - d1) != 0)
+		sol -= (x1 - y1 - d1) * g0(0.5, -y1 - d1, -x1,
+				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));
+
+	if ((x2 - y2 - d2) != 0)
+		sol -= (x2 - y2 - d2) * g0(0.5, x2, y2 + d2,
+				sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));
+
+	double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2
+			- d2) + d3 * d3);
+	sol += 1. / 3 * hlp * sqrt(hlp);
+	return sol;
+}
+
 double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,
 		double d3) {
 
diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp
index 26ebd71..87e3bdf 100644
--- a/src/slpRectangle.hpp
+++ b/src/slpRectangle.hpp
@@ -9,11 +9,20 @@ double g0(double, double, double, double);
 // sol = G00(p,y1,y2,x1,x2,l);
 double G00(double, double, double, double, double, double);
 // sol = F_par(x1,x2,y1,y2,d1,d2,d3);
-double F_par(double, double, double, double, double, double, double);
+//double F_par(double, double, double, double, double, double, double);
+#define F_par F_parY1Y2
 // sol = F_ort(x1,x2,y2,y3,d1,d2,d3);
 double F_ort(double, double, double, double, double, double, double);
 
 
+
+// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
+double F_parY1Y2(double, double, double, double, double, double, double);
+// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
+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);
+
 double slpADLO(double, double, double, double, double);
 double compute_g0(double, double, double, double);
 double FLO_plane(double, double, double, double, double, double, double);
diff --git a/src/test_solveError.m b/src/test_solveError.m
index b9a94a5..0316541 100644
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -2,8 +2,8 @@
 dataIso =[];
 dataAniso =[];
 maxSize = 1050;
-% Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
-Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren
+Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
+% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren
 % Netz = 'exmpl_2DQuad';
 
 %% Isotrope Verfeinerung
@@ -19,12 +19,12 @@ while size(elements,1)<maxSize
   A_fine = mex_build_A(coordinates_fine,elements_fine);
   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
   x_fine = A_fine\b;
-%   xe_fine = x_fine'*A_fine*x_fine;
+  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))];
+  dataIso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
   [coordinates, elements] = refineQuad(coordinates,elements,marked);
   toc
 end
@@ -44,12 +44,12 @@ while size(elements,1)<maxSize
   A_fine = mex_build_A(coordinates_fine,elements_fine);
   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
   x_fine = A_fine\b;
-%   xe_fine = x_fine'*A_fine*x_fine;
+  xe_fine = x_fine'*A_fine*x_fine;
 
   ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
   marked = mark(x_fine(f2s)',ind,0.5,0.5);
 
-  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind))];
+  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
   
   if(sum(marked~=1)==0)
     disp 'Abgebrochen weil nicht mehr verfeinert wird!'
@@ -65,8 +65,16 @@ 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