[src] Analytisch Y1X2 eingebaut sowie X1Y2 -> für für paralelle Elemente

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

diff --git a/src/mex_FparY1X2.cpp b/src/mex_FparY1X2.cpp
new file mode 100755 (executable)
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 7e9a7b0a6894192d4effd1c83369b02b1a0b90b3..42f480a3b56d28efc599fb21f770b45a91c9a401 100644 (file)
--- 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;\r
 }\r
 \r
-double F_par(double x1, double x2, double y1, double y2, double d1, double d2,\r
+double F_parY1Y2(double x1, double x2, double y1, double y2, double d1, double d2,\r
                double d3) {\r
 \r
        //       printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
@@ -109,6 +109,52 @@ double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
        return sol;\r
 }\r
 \r
+double F_parY1X2(double x1, double x2, double y1, double y2, double d1, double d2,\r
+               double d3) {\r
+\r
+       //       printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
+       double sol = (x1 - y1 - d1) * (x2 - y2 - d2);\r
+\r
+       if (sol != 0)\r
+               sol *= G00(-0.5, x1, -y2 - d2, y1 + d1, -x2, d3);\r
+\r
+       if ((x1 - y1 - d1) != 0)\r
+               sol -= (x1 - y1 - d1) * g0(0.5, x1, y1 + d1,\r
+                               sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));\r
+\r
+       if ((x2 - y2 - d2) != 0)\r
+               sol -= (x2 - y2 - d2) * g0(0.5, -y2 - d2, -x2,\r
+                               sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));\r
+\r
+       double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
+                       - d2) + d3 * d3);\r
+       sol += 1. / 3 * hlp * sqrt(hlp);\r
+       return sol;\r
+}\r
+\r
+double F_parX1Y2(double x1, double x2, double y1, double y2, double d1, double d2,\r
+               double d3) {\r
+\r
+       //       printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
+       double sol = (x1 - y1 - d1) * (x2 - y2 - d2);\r
+\r
+       if (sol != 0)\r
+               sol *= G00(-0.5, -y1 - d1, x2, -x1, y2 + d2, d3);\r
+\r
+       if ((x1 - y1 - d1) != 0)\r
+               sol -= (x1 - y1 - d1) * g0(0.5, -y1 - d1, -x1,\r
+                               sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));\r
+\r
+       if ((x2 - y2 - d2) != 0)\r
+               sol -= (x2 - y2 - d2) * g0(0.5, x2, y2 + d2,\r
+                               sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));\r
+\r
+       double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
+                       - d2) + d3 * d3);\r
+       sol += 1. / 3 * hlp * sqrt(hlp);\r
+       return sol;\r
+}\r
+\r
 double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,\r
                double d3) {\r
 \r

diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp
index 26ebd71a3da57e63b08423cc5f75f3bbf653eb36..87e3bdf24b11a6b7e95bba3c7c4420d98bafc280 100644 (file)
--- 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 b9a94a55f5a8b648503b180c3eeb4e3d89a62a67..0316541b8239c87108c54fd3fc525a2632731bad 100644 (file)
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -2,8 +2,8 @@
 dataIso =[];\r
 dataAniso =[];\r
 maxSize = 1050;\r
-% Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren\r
-Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren\r
+Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren\r
+% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren\r
 % Netz = 'exmpl_2DQuad';\r
 \r
 %% Isotrope Verfeinerung\r
@@ -19,12 +19,12 @@ while size(elements,1)<maxSize
   A_fine = mex_build_A(coordinates_fine,elements_fine);\r
   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
   x_fine = A_fine\b;\r
-%   xe_fine = x_fine'*A_fine*x_fine;\r
+  xe_fine = x_fine'*A_fine*x_fine;\r
 \r
   ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
   marked = 2*ones(1,size(elements,1));\r
 \r
-  dataIso(ref,:) = [size(elements,1) sqrt(sum(ind))];\r
+  dataIso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
   [coordinates, elements] = refineQuad(coordinates,elements,marked);\r
   toc\r
 end\r
@@ -44,12 +44,12 @@ while size(elements,1)<maxSize
   A_fine = mex_build_A(coordinates_fine,elements_fine);\r
   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
   x_fine = A_fine\b;\r
-%   xe_fine = x_fine'*A_fine*x_fine;\r
+  xe_fine = x_fine'*A_fine*x_fine;\r
 \r
   ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
   marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
 \r
-  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind))];\r
+  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
   \r
   if(sum(marked~=1)==0)\r
     disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
@@ -65,8 +65,16 @@ clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine
 \r
 \r
 %% Alles Zeichnen\r
+figure (1)\r
 loglog(dataIso(:,1),dataIso(:,2),dataAniso(:,1),dataAniso(:,2));\r
 \r
 legend('Isotrop','Anisotrop')\r
 xlabel('Elemente')\r
 ylabel('Fehler')\r
+\r
+figure (2)\r
+loglog(dataIso(:,1),dataIso(:,3),dataAniso(:,1),dataAniso(:,3));\r
+\r
+legend('Isotrop','Anisotrop')\r
+xlabel('Elemente')\r
+ylabel('EnergieNorm')\r