[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

diff --git a/src/mex_Fpar.cpp b/src/mex_Fpar.cpp
index c02cfb60d095319e29f2c7cd8fcc400e6625a7ed..9103e4d6f83c9228963ec66c9c95a1d417d7e70b 100644 (file)
--- 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 9f9fd9876c28c6087f378aca9f65fd7386e1b126..24379f508f916f4a79bf37ef2c0404fa75ec1943 100644 (file)
--- 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 48a30b0de4f8209b260093bd2ad41baba60cbfd0..0cdec5f06a67e50434c13b5db376050493c3aae0 100644 (file)
--- a/src/mex_build_A.cpp
+++ b/src/mex_build_A.cpp
@@ -1,3 +1,11 @@
+/*\r
+ * voll Analytisch\r
+ * ohne Optimierungen\r
+ *\r
+ *\r
+ */\r
+\r
+\r
 #include <iostream>\r
 #include <cmath>\r
 #include <cassert>\r
@@ -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]);\r
                                        //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
                                        tmp\r
-                                                       = apply0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),\r
+                                                       = apply0Int4(F_par, fabs(xa[rxa]), fabs(xb[rxb]),\r
                                                                        fabs(ya[rxa]), fabs(yb[rxb]), d[rxa],\r
                                                                        d[rxb], d[rx]);\r
                                        //       printf("%d%d|%.2f\n",j,k,tmp);\r
@@ -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]);\r
                                        //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
                                        tmp\r
-                                                       = apply0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),\r
+                                                       = apply0Int4(F_par, fabs(xa[rxa]), fabs(xb[rxb]),\r
                                                                        fabs(yb[rxa]), fabs(ya[rxb]), d[rxa],\r
                                                                        d[rxb], d[rx]);\r
                                        //       printf("%d%d|%.2f\n",j,k,tmp);\r
@@ -247,22 +255,22 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 \r
                                if (rxa == rya) {\r
                                        tmp\r
-                                                       = apply0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),\r
+                                                       = apply0Int4(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),\r
                                                                        fabs(ya[rya]), fabs(yb[ryb]), d[rxb],\r
                                                                        d[rxa], d[rx]);\r
                                } else if (rxa == ryb) {\r
                                        tmp\r
-                                                       = apply0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),\r
+                                                       = apply0Int4(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),\r
                                                                        fabs(yb[ryb]), fabs(ya[rya]), d[rxb],\r
                                                                        d[rxa], d[rx]);\r
                                } else if (rxb == rya) {\r
                                        tmp\r
-                                                       = apply0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),\r
+                                                       = apply0Int4(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),\r
                                                                        fabs(ya[rya]), fabs(yb[ryb]), d[rxa],\r
                                                                        d[rxb], d[rx]);\r
                                } else {\r
                                        tmp\r
-                                                       = apply0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),\r
+                                                       = apply0Int4(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),\r
                                                                        fabs(yb[ryb]), fabs(ya[rya]), d[rxa],\r
                                                                        d[rxb], d[rx]);\r
                                }\r

diff --git a/src/mex_build_As1.cpp b/src/mex_build_As1.cpp
index 485b74fb199f9a2ac2b4ac337c6412896cd8b5c0..2356c72c3b12cd5f1f6dc2078a199c9931289a16 100644 (file)
--- a/src/mex_build_As1.cpp
+++ b/src/mex_build_As1.cpp
@@ -1,3 +1,11 @@
+/*\r
+ * voll Analytisch\r
+ * Element mit groesserer Diagonale nach aussen gedreht\r
+ *\r
+ *\r
+ */\r
+\r
+\r
 #include <iostream>\r
 #include <cmath>\r
 #include <cassert>\r
@@ -72,11 +80,14 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
        double tmp,tmp2;\r
 \r
        int itmp;\r
-\r
+       double dtmp;\r
        double * vtmp;\r
 \r
+       long double xda,xdb,yda,ydb;\r
+\r
        //LageInformationen\r
        int rx, rxa, rxb, ry, rya, ryb;\r
+       int rtxa, rtxb, rtx;\r
 \r
        //Ausrechnen\r
        for (j = 0; j < em; ++j) {\r
@@ -110,17 +121,17 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                // Lageeigenschaften des Flaechenstuecks\r
                //if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])\r
 \r
-               rx = dimOfVec(xn);\r
+               rtx = dimOfVec(xn);\r
 \r
-               rxa = dimOfVec(xa);\r
+               rtxa = dimOfVec(xa);\r
 \r
-               rxb = dimOfVec(xb);\r
+               rtxb = dimOfVec(xb);\r
 \r
 \r
                //kleinste Ecke finden und fuer \delta verwenden\r
 \r
-               if (xa[rxa] > 0) {\r
-                       if (xb[rxb] > 0) {\r
+               if (xa[rtxa] > 0) {\r
+                       if (xb[rtxb] > 0) {\r
                                for (i = 0; i < 3; ++i)\r
                                        dt[i] = -x0[i];\r
                        } else {\r
@@ -128,7 +139,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                                        dt[i] = -x3[i];\r
                        }\r
                } else {\r
-                       if (xb[rxb] > 0) {\r
+                       if (xb[rtxb] > 0) {\r
                                for (i = 0; i < 3; ++i)\r
                                        dt[i] = -x1[i];\r
                        } else {\r
@@ -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]);\r
 \r
                for (k = 0; k < em; ++k) {\r
+\r
+                       xda = xa[rtxa];\r
+                       xdb = xb[rtxb];\r
+\r
+                       rx = rtx;\r
+                       rxa = rtxa;\r
+                       rxb = rtxb;\r
+\r
+\r
                        y0[0] = C[(int) E[k] - 1];\r
                        y0[1] = C[cm + (int) E[k] - 1];\r
                        y0[2] = C[2 * cm + (int) E[k] - 1];\r
@@ -198,22 +218,20 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                                }\r
                        }\r
 \r
-                       //      for (i = 0; i<3;++i)\r
-                       //              d[i] = y0[i]-x0[i];\r
-\r
+                       yda = ya[rya];\r
+                       ydb = yb[ryb];\r
 \r
-                       if(ry == rx){\r
 \r
-                               if(xa[rxa]*xa[rxa]+xb[rxb]*xb[rxb]>ya[rya]*ya[rya]+yb[ryb]*yb[ryb]){\r
+                               if(xda*xda+xdb*xdb>yda*yda+ydb*ydb){\r
 \r
-                                       printf("+(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
-                                       printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);\r
+                               /*      printf("+(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
+                                       printf("[%.2f,%.2f@%.2f,%.2f]",xda,xb[rxb],ya[rya],yb[ryb]);\r
                                        printf("{%d%d} %.1f %.1f %.1f",j,k,d[0],d[1],d[2]);\r
-                                       printf("\n");\r
+                                       printf("\n");*/\r
 \r
 \r
-                                       vtmp = xa; xa = ya; ya = vtmp;\r
-                                       vtmp = xb; xb = yb; yb = vtmp;\r
+                                       dtmp = xda; xda = yda; yda = dtmp;\r
+                                       dtmp = xdb; xdb = ydb; ydb = dtmp;\r
 \r
                                        itmp = rxa; rxa = rya; rya = itmp;\r
                                        itmp = rxb; rxb = ryb; ryb = itmp;\r
@@ -224,19 +242,17 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 \r
 \r
 \r
-                                       printf("-(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
-                                       printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);\r
+                       /*              printf("-(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
+                                       printf("[%.2f,%.2f@%.2f,%.2f]",xda,xb[rxb],ya[rya],yb[ryb]);\r
                                        printf("{%d%d} %.1f %.1f %.1f",j,k,d[0],d[1],d[2]);\r
-                                       printf("\n");\r
+                                       printf("\n");*/\r
                                }\r
 \r
 \r
-                       }\r
-\r
                        //printf("(%d,%d)",rx,ry);\r
                        if (rx == ry) {\r
                                //printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
-                               //printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);\r
+                               //printf("[%.2f,%.2f@%.2f,%.2f]",xda,xdb,yda,ydb);\r
                                //printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);\r
                                //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
                                if (rxa == rya) {\r
@@ -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]);\r
                                        //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
 \r
-                                       if(xa[rxa]*xa[rxa]+xb[rxb]*xb[rxb]<=ya[rya]*ya[rya]+yb[ryb]*yb[ryb]){\r
+                       //              if(xda*xda+xdb*xdb<=yda*yda+ydb*ydb){\r
 \r
                                        tmp\r
-                                                       = apply0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),\r
-                                                                       fabs(ya[rya]), fabs(yb[ryb]), d[rxa],\r
+                                                       = apply0Int4(F_par, fabs(xda), fabs(xdb),\r
+                                                                       fabs(yda), fabs(ydb), d[rxa],\r
                                                                        d[rxb], d[rx]);\r
 \r
-                                       }else{\r
+                       //              }else{\r
 \r
-                                       tmp\r
-                                                       = apply0Int(F_par, fabs(ya[rya]), fabs(yb[ryb]),\r
-                                                                       fabs(xa[rxa]), fabs(xb[rxb]), -d[rxa],\r
-                                                                       -d[rxb], -d[rx]);\r
-                                       }\r
+                       //              tmp\r
+                       //                              = apply0Int4(F_par, fabs(yda), fabs(ydb),\r
+                       //                                              fabs(xda), fabs(xdb), -d[rxa],\r
+                       //                                              -d[rxb], -d[rx]);\r
+                       //              }\r
                        //              if(fabs(tmp-tmp2)>10e-16)\r
                        //                      printf("wtf");\r
                                        //       printf("%d%d|%.2f\n",j,k,tmp);\r
@@ -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]);\r
                                        //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
                                        tmp\r
-                                                       = apply0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),\r
-                                                                       fabs(yb[rxa]), fabs(ya[rxb]), d[rxa],\r
+                                                       = apply0Int4(F_par, fabs(xda), fabs(xdb),\r
+                                                                       fabs(ydb), fabs(yda), d[rxa],\r
                                                                        d[rxb], d[rx]);\r
                                        //       printf("%d%d|%.2f\n",j,k,tmp);\r
                                }\r
@@ -279,26 +295,25 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 \r
                                if (rxa == rya) {\r
                                        tmp\r
-                                                       = apply0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),\r
-                                                                       fabs(ya[rya]), fabs(yb[ryb]), d[rxb],\r
+                                                       = apply0Int4(F_ort, fabs(xdb), fabs(xda),\r
+                                                                       fabs(yda), fabs(ydb), d[rxb],\r
                                                                        d[rxa], d[rx]);\r
                                } else if (rxa == ryb) {\r
                                        tmp\r
-                                                       = apply0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),\r
-                                                                       fabs(yb[ryb]), fabs(ya[rya]), d[rxb],\r
+                                                       = apply0Int4(F_ort, fabs(xdb), fabs(xda),\r
+                                                                       fabs(ydb), fabs(yda), d[rxb],\r
                                                                        d[rxa], d[rx]);\r
                                } else if (rxb == rya) {\r
                                        tmp\r
-                                                       = apply0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),\r
-                                                                       fabs(ya[rya]), fabs(yb[ryb]), d[rxa],\r
+                                                       = apply0Int4(F_ort, fabs(xda), fabs(xdb),\r
+                                                                       fabs(yda), fabs(ydb), d[rxa],\r
                                                                        d[rxb], d[rx]);\r
                                } else {\r
                                        tmp\r
-                                                       = apply0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),\r
-                                                                       fabs(yb[ryb]), fabs(ya[rya]), d[rxa],\r
+                                                       = apply0Int4(F_ort, fabs(xda), fabs(xdb),\r
+                                                                       fabs(ydb), fabs(yda), d[rxa],\r
                                                                        d[rxb], d[rx]);\r
                                }\r
-                               //A[(k * em) + j] = tmp;\r
 \r
                        }\r
                        A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r

diff --git a/src/mex_build_As2.cpp b/src/mex_build_As2.cpp
new file mode 100644 (file)
index 0000000..c602509
--- /dev/null
+++ b/src/mex_build_As2.cpp
@@ -0,0 +1,328 @@
+/*\r
+ * voll Analytisch\r
+ * Element mit groesserer Diagonale nach aussen gedreht\r
+ *\r
+ *\r
+ */\r
+\r
+\r
+#include <iostream>\r
+#include <cmath>\r
+#include <cassert>\r
+#include "mex.h"\r
+#include <stdlib.h>\r
+\r
+//#include "tbb/parallel_for.h"\r
+\r
+\r
+#include "slpRectangle.hpp"\r
+\r
+#define EPS 0.02\r
+\r
+using namespace std;\r
+\r
+inline int dimOfVec(double * vec)\r
+{\r
+       if (vec[2] != 0)\r
+               return 2;\r
+       else if (vec[1] != 0)\r
+               return 1;\r
+       else\r
+               return 0;\r
+}\r
+\r
+\r
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
+\r
+       int i, j, k;\r
+\r
+       //sicherheitsabfragen zu Datengroessen\r
+       if (nrhs != 2)\r
+               mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx4))");\r
+       if (nlhs > 1)\r
+               mexErrMsgTxt("has only one output argument");\r
+\r
+       int cm = mxGetM(prhs[0]);\r
+       int cn = mxGetN(prhs[0]);\r
+\r
+       if (cn != 3)\r
+               mexErrMsgTxt("expected coordinates (Nx3)");\r
+       int em = mxGetM(prhs[1]);\r
+       int en = mxGetN(prhs[1]);\r
+       if (en != 4)\r
+               mexErrMsgTxt("expected elements (Mx4)");\r
+       //Vorbereitung der Daten\r
+\r
+       plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);\r
+       double * A = mxGetPr(plhs[0]);\r
+       double * C = mxGetPr(prhs[0]);\r
+       double * E = mxGetPr(prhs[1]);\r
+\r
+       double * x0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * xn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+\r
+       double * y0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * yn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+\r
+       double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+       double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+\r
+       double tmp,tmp2;\r
+\r
+       int itmp;\r
+       double dtmp;\r
+       double * vtmp;\r
+\r
+       long double xda,xdb,yda,ydb;\r
+\r
+       //LageInformationen\r
+       int rx, rxa, rxb, ry, rya, ryb;\r
+       int rtxa, rtxb, rtx;\r
+\r
+       //Ausrechnen\r
+       for (j = 0; j < em; ++j) {\r
+               x0[0] = C[(int) E[j] - 1];\r
+               x0[1] = C[cm + (int) E[j] - 1];\r
+               x0[2] = C[2 * cm + (int) E[j] - 1];\r
+\r
+               x1[0] = C[(int) E[em + j] - 1];\r
+               x1[1] = C[cm + (int) E[em + j] - 1];\r
+               x1[2] = C[2 * cm + (int) E[em + j] - 1];\r
+\r
+               x2[0] = C[(int) E[2 * em + j] - 1];\r
+               x2[1] = C[cm + (int) E[2 * em + j] - 1];\r
+               x2[2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
+\r
+               x3[0] = C[(int) E[3 * em + j] - 1];\r
+               x3[1] = C[cm + (int) E[3 * em + j] - 1];\r
+               x3[2] = C[2 * cm + (int) E[3 * em + j] - 1];\r
+\r
+               for (i = 0; i < 3; ++i)\r
+                       xa[i] = x1[i] - x0[i];\r
+\r
+               for (i = 0; i < 3; ++i)\r
+                       xb[i] = x3[i] - x0[i];\r
+\r
+               //Kreuzprodukt\r
+               xn[0] = (xa[1]) * (xb[2]) - (xa[2]) * (xb[1]);\r
+               xn[1] = (xa[2]) * (xb[0]) - (xa[0]) * (xb[2]);\r
+               xn[2] = (xa[0]) * (xb[1]) - (xa[1]) * (xb[0]);\r
+\r
+               // Lageeigenschaften des Flaechenstuecks\r
+               //if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])\r
+\r
+               rtx = dimOfVec(xn);\r
+\r
+               rtxa = dimOfVec(xa);\r
+\r
+               rtxb = dimOfVec(xb);\r
+\r
+\r
+               //kleinste Ecke finden und fuer \delta verwenden\r
+\r
+               if (xa[rtxa] > 0) {\r
+                       if (xb[rtxb] > 0) {\r
+                               for (i = 0; i < 3; ++i)\r
+                                       dt[i] = -x0[i];\r
+                       } else {\r
+                               for (i = 0; i < 3; ++i)\r
+                                       dt[i] = -x3[i];\r
+                       }\r
+               } else {\r
+                       if (xb[rtxb] > 0) {\r
+                               for (i = 0; i < 3; ++i)\r
+                                       dt[i] = -x1[i];\r
+                       } else {\r
+                               for (i = 0; i < 3; ++i)\r
+                                       dt[i] = -x2[i];\r
+                               //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
+                       }\r
+               }\r
+\r
+               //for (i=0;i<3;++i)\r
+               //      dt[i] = 0;\r
+\r
+\r
+               // printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);\r
+\r
+               for (k = 0; k < em; ++k) {\r
+\r
+                       xda = xa[rtxa];\r
+                       xdb = xb[rtxb];\r
+\r
+                       rx = rtx;\r
+                       rxa = rtxa;\r
+                       rxb = rtxb;\r
+\r
+\r
+                       y0[0] = C[(int) E[k] - 1];\r
+                       y0[1] = C[cm + (int) E[k] - 1];\r
+                       y0[2] = C[2 * cm + (int) E[k] - 1];\r
+\r
+                       y1[0] = C[(int) E[em + k] - 1];\r
+                       y1[1] = C[cm + (int) E[em + k] - 1];\r
+                       y1[2] = C[2 * cm + (int) E[em + k] - 1];\r
+\r
+                       y2[0] = C[(int) E[2 * em + k] - 1];\r
+                       y2[1] = C[cm + (int) E[2 * em + k] - 1];\r
+                       y2[2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
+\r
+                       y3[0] = C[(int) E[3 * em + k] - 1];\r
+                       y3[1] = C[cm + (int) E[3 * em + k] - 1];\r
+                       y3[2] = C[2 * cm + (int) E[3 * em + k] - 1];\r
+\r
+                       for (i = 0; i < 3; ++i)\r
+                               ya[i] = y1[i] - y0[i];\r
+\r
+                       for (i = 0; i < 3; ++i)\r
+                               yb[i] = y3[i] - y0[i];\r
+\r
+                       yn[0] = (ya[1]) * (yb[2]) - (ya[2]) * (yb[1]);\r
+                       yn[1] = (ya[2]) * (yb[0]) - (ya[0]) * (yb[2]);\r
+                       yn[2] = (ya[0]) * (yb[1]) - (ya[1]) * (yb[0]);\r
+\r
+                       ry = dimOfVec(yn);\r
+\r
+                       rya = dimOfVec(ya);\r
+\r
+                       ryb = dimOfVec(yb);\r
+\r
+                       //kleinste Ecke finden und fuer \delta verwenden\r
+\r
+                       if (ya[rya] > 0) {\r
+                               if (yb[ryb] > 0) {\r
+                                       for (i = 0; i < 3; ++i)\r
+                                               d[i] = dt[i] + y0[i];\r
+                               } else {\r
+                                       for (i = 0; i < 3; ++i)\r
+                                               d[i] = dt[i] + y3[i];\r
+                               }\r
+                       } else {\r
+                               if (yb[ryb] > 0) {\r
+                                       for (i = 0; i < 3; ++i)\r
+                                               d[i] = dt[i] + y1[i];\r
+                               } else {\r
+                                       for (i = 0; i < 3; ++i)\r
+                                               d[i] = dt[i] + y2[i];\r
+                                       //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
+                               }\r
+                       }\r
+\r
+                       yda = ya[rya];\r
+                       ydb = yb[ryb];\r
+\r
+                               //groesseres Element nach aussen drehen\r
+                               if(xda*xda+xdb*xdb>yda*yda+ydb*ydb){\r
+\r
+                                       dtmp = xda; xda = yda; yda = dtmp;\r
+                                       dtmp = xdb; xdb = ydb; ydb = dtmp;\r
+\r
+                                       itmp = rxa; rxa = rya; rya = itmp;\r
+                                       itmp = rxb; rxb = ryb; ryb = itmp;\r
+                                       itmp = rx; rx = ry; ry = itmp;\r
+\r
+                                       for (i = 0; i<3;++i)\r
+                                               d[i] = -d[i];\r
+\r
+                               }\r
+\r
+                               //if(xda*xda+xdb*xdb>=d[0]*d[0]+d[1]*d[1]+d[2]*d[2]){\r
+                               //      printf(".");\r
+                               //}\r
+\r
+                       //printf("(%d,%d)",rx,ry);\r
+                       if (rx == ry) {\r
+                               //printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
+                               //printf("[%.2f,%.2f@%.2f,%.2f]",xda,xdb,yda,ydb);\r
+                               //printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);\r
+                               //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
+                               if (rxa == rya) {\r
+                                       //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],\r
+                                       //              y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
+                                       //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
+\r
+                       //              if(xda*xda+xdb*xdb<d[0]*d[0]+d[1]*d[1]+d[2]*d[2]){\r
+\r
+                                       tmp\r
+                                                       = apply0Int4(F_par, fabs(xda), fabs(xdb),\r
+                                                                       fabs(yda), fabs(ydb), d[rxa],\r
+                                                                       d[rxb], d[rx]);\r
+\r
+                       //              }else{\r
+\r
+                                       tmp2\r
+                                                       = apply0Int2(F_parY1Y2_2, fabs(xda), fabs(xdb),\r
+                                                                       fabs(yda), fabs(ydb), d[rxa],\r
+                                                                       d[rxb], d[rx]);\r
+                       //              }\r
+                                       if(fabs(tmp-tmp2)>10e-8)\r
+                                               return;\r
+                                       //       printf("%d%d|%.2f\n",j,k,tmp);\r
+                               } else {\r
+                                       //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],\r
+                                       //              y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
+                                       //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
+                                       tmp\r
+                                                       = apply0Int4(F_par, fabs(xda), fabs(xdb),\r
+                                                                       fabs(ydb), fabs(yda), d[rxa],\r
+                                                                       d[rxb], d[rx]);\r
+                                       //       printf("%d%d|%.2f\n",j,k,tmp);\r
+                               }\r
+\r
+                       } else {\r
+                               //printf("(%d,%d@%d,%d)", rxa, rxb, rya, ryb);\r
+                               //printf("(%d,%d)", rx, ry);\r
+                               //printf("\n");\r
+\r
+                               if (rxa == rya) {\r
+                                       tmp\r
+                                                       = apply0Int4(F_ort, fabs(xdb), fabs(xda),\r
+                                                                       fabs(yda), fabs(ydb), d[rxb],\r
+                                                                       d[rxa], d[rx]);\r
+                               } else if (rxa == ryb) {\r
+                                       tmp\r
+                                                       = apply0Int4(F_ort, fabs(xdb), fabs(xda),\r
+                                                                       fabs(ydb), fabs(yda), d[rxb],\r
+                                                                       d[rxa], d[rx]);\r
+                               } else if (rxb == rya) {\r
+                                       tmp\r
+                                                       = apply0Int4(F_ort, fabs(xda), fabs(xdb),\r
+                                                                       fabs(yda), fabs(ydb), d[rxa],\r
+                                                                       d[rxb], d[rx]);\r
+                               } else {\r
+                                       tmp\r
+                                                       = apply0Int4(F_ort, fabs(xda), fabs(xdb),\r
+                                                                       fabs(ydb), fabs(yda), d[rxa],\r
+                                                                       d[rxb], d[rx]);\r
+                               }\r
+\r
+                       }\r
+                       A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
+\r
+               }\r
+\r
+       }\r
+       //printf("\n");\r
+       //Rueckgabe (eventuell zurueck schreiben)\r
+       //   mxFree(x0);\r
+       //mxFree(x1);\r
+       //mxFree(x3);\r
+       //mxFree(y0);\r
+       //mxFree(y1);\r
+       //mxFree(y3);\r
+       //mxFree(xn);\r
+       //mxFree(yn);\r
+       //mxFree(d);\r
+\r
+       return;\r
+}\r

diff --git a/src/mex_g0.cpp b/src/mex_g0.cpp
index 992481361ee6c7c811808c1126959882efc7682c..7dfb788555cf72c583bca28c9be2bc705898b893 100644 (file)
--- 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 ed3eaad3762526d825fd0f16f18ae2148da7c2bd..f868508cb42165a61125148f8ae48db23eb6b7c1 100644 (file)
--- 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 42f480a3b56d28efc599fb21f770b45a91c9a401..b92c6a40093524eaa3f984d27ea15ae4c28ff016 100644 (file)
--- a/src/slpRectangle.cpp
+++ b/src/slpRectangle.cpp
@@ -3,12 +3,17 @@
 #include <cassert>\r
 #include <stdlib.h>\r
 \r
+#include <mex.h>\r
+\r
 #include "slpRectangle.hpp"\r
 \r
-#define EPS 0.00001\r
+#define EPS 10e5\r
 \r
 using namespace std;\r
 \r
+double gauss_w5[] = {0.1185,0.2393,0.2844,0.2393,0.1185};\r
+double gauss_n5[] = {0.0469,0.2308,0.5000,0.7692,0.9531};\r
+\r
 int sign(double);\r
 //double arsinh(double);\r
 \r
@@ -20,6 +25,31 @@ int inline sign(double x) {
  return log(x + sqrt(x * x + 1));\r
  }\r
  */\r
+\r
+double singleQuad(double p, double y, double x, double l) {\r
+\r
+       double sol = 0;\r
+\r
+       for(int i = 0;i<5;++i){\r
+               sol += pow((x-gauss_n5[i]*y)*(x-gauss_n5[i]*y)+l*l,p)*gauss_w5[i]*y;\r
+       }\r
+\r
+       return sol;\r
+}\r
+\r
+\r
+double doubleQuad(double p, double y1, double y2, double x1, double x2, double l) {\r
+\r
+       double sol = 0;\r
+       for(int i = 0;i<5;++i){\r
+               for(int j=0;j<5;++j){\r
+                       sol += pow((x1-y1*gauss_n5[i])*(x1-y1*gauss_n5[i])\r
+                                       +(x2-y2*gauss_n5[j])*(x2-y2*gauss_n5[j])+l*l,p)*y1*gauss_w5[i]*y2*gauss_w5[j];\r
+               }\r
+       }\r
+       return sol;\r
+}\r
+\r
 //y-x muss != 0 sein\r
 double g0(double p, double y, double x, double l) {\r
        //printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);\r
@@ -155,6 +185,29 @@ double F_parX1Y2(double x1, double x2, double y1, double y2, double d1, double d
        return sol;\r
 }\r
 \r
+double F_parY1Y2_2(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 *= doubleQuad(-0.5, x1, x2, y1 + d1, y2 + d2, d3);\r
+\r
+       if ((x1 - y1 - d1) != 0)\r
+               sol -= (x1 - y1 - d1) * singleQuad(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) * singleQuad(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
@@ -180,7 +233,7 @@ double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,
        return sol / 2.;\r
 }\r
 \r
-double applyInt(\r
+double applyInt4(\r
                double(*f)(double, double, double, double, double, double, double),\r
                double a, double b, double c, double d, double r, double t, double u,\r
                double v, double d1, double d2, double d3) {\r
@@ -205,7 +258,7 @@ double applyInt(
         */\r
 }\r
 \r
-double apply0Int(\r
+double apply0Int4(\r
                double(*f)(double, double, double, double, double, double, double),\r
                double b, double d, double t, double v, double d1, double d2, double d3) {\r
 \r
@@ -219,6 +272,17 @@ double apply0Int(
 \r
 }\r
 \r
+double apply0Int2(\r
+               double(*f)(double, double, double, double, double, double, double),\r
+               double b, double d, double t, double v, double d1, double d2, double d3) {\r
+\r
+       return f(b, d, t, v, d1, d2, d3)\r
+                       -f(b, 0, t, v, d1, d2, d3)\r
+                       -f(0, d, t, v, d1, d2, d3)\r
+                       +f(0, 0, t, v, d1, d2, d3);\r
+\r
+}\r
+\r
 double slpADLO(double y1, double y2, double x1, double x2, double a) {\r
        double G3 = 0;\r
        double gL = 0;\r

diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp
index 87e3bdf24b11a6b7e95bba3c7c4420d98bafc280..5b3a715ec35acf6f12ba7b65ac6a9ad6326c8b64 100644 (file)
--- 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 27c3e6c559898d68312a26755cd0f169f798d862..e1f70a7ae7575105dfbca182c3832852ed1498f9 100644 (file)
--- 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 c23354f013e223601a020f3033dda1611a1b3f9f..66e7f7af4a382be85e4a3ccd5691010951c0c990 100644 (file)
--- 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 89e2a02a8f80a0078b5caa444cd60ab512c37c32..dbd75dcfa00c73969f87edac4842109a2c282989 100644 (file)
--- 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 c367bb2e107fc67b71df613a97ee4648deeca192..b165322a6743046b939d4321848666d7830efd85 100755 (executable)
--- a/src/test_solveErrorS1.m
+++ b/src/test_solveErrorS1.m
@@ -3,7 +3,7 @@ mex mex_build_As1.cpp slpRectangle.cpp
 \r
 dataIso =[];\r
 dataAniso =[];\r
-maxSize = 10^2;\r
+maxSize = 10^3;\r
 Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren\r
 % Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren\r
 % Netz = 'exmpl_2DQuad';\r
@@ -47,23 +47,23 @@ while size(elements,1)<maxSize
   \r
   %Matrix mit feinem Netz berechnen\r
   A_fine = mex_build_As1(coordinates_fine,elements_fine);\r
-  A2 = mex_build_A(coordinates_fine,elements_fine);\r
+%   A2 = mex_build_A(coordinates_fine,elements_fine);\r
   \r
 %   figure(1)\r
 %   spy(A_fine-A_fine')\r
 %   \r
 %   figure (2)\r
 %   spy(A2-A2')\r
-  eh = A_fine - A2;\r
+%   eh = A_fine - A2;\r
 %   sum(sum(abs(A2-A2')))\r
 %   sum(sum(abs(A_fine-A_fine')))\r
   \r
-  spy(eh)\r
+%   spy(eh)\r
  \r
-  sum(sum(abs(eh)))\r
-  if(eh.^2)\r
-    disp 'hä'\r
-  end\r
+%   sum(sum(abs(eh)))\r
+%   if(eh.^2)\r
+%     disp 'hä'\r
+%   end\r
   \r
   %rechte Seite Aufbauen\r
   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r

diff --git a/src/test_solveError_CGrid.m b/src/test_solveError_CGrid.m
new file mode 100755 (executable)
index 0000000..4c745b6
--- /dev/null
+++ b/src/test_solveError_CGrid.m
@@ -0,0 +1,57 @@
+\r
+dataIso =[];\r
+dataAniso =[];\r
+maxSize = 500;\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
+%% Anisotrope Verfeinerung\r
+disp 'Anisotrop'\r
+load(Netz)\r
+\r
+ref = 0;\r
+while size(elements,1)<maxSize\r
+  tic\r
+  ref = ref+1; \r
+\r
+  %Netz Isotrop verfeinern (4 Teile)\r
+  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
+  \r
+  %Matrix mit feinem Netz berechnen\r
+  A_fine = mex_build_A(coordinates_fine,elements_fine);\r
+  \r
+  %rechte Seite Aufbauen\r
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
+  \r
+  %Lösungsvektor bestimmen\r
+  x_fine = A_fine\b;\r
+  \r
+  %Energienorm^2 bestimmen\r
+  xe_fine = x_fine'*A_fine*x_fine;\r
+\r
+  % Welche Elemente sollen verfeinert werden\r
+  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
+  \r
+  % Wie soll verfeinert werden\r
+  marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
+\r
+  %Daten zur Auswertung zwischenspeichern (testweise?)\r
+  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
+  \r
+  %Sicherheit falls keine Elemente markiert sind\r
+  if(sum(marked~=1)==0)\r
+    disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
+    break;\r
+  end\r
+  \r
+  %Verfeinern\r
+  [coordinates, elements] = refineQuad(coordinates,elements,marked);\r
+  file = ['save' int2str(ref)];\r
+  save ([file], 'coordinates', 'elements')\r
+  toc\r
+end\r
+dataAniso\r
+\r
+clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine\r
+\r

diff --git a/src/test_solveError_Grid.m b/src/test_solveError_Grid.m
deleted file mode 100755 (executable)
index 4c745b6..0000000
--- a/src/test_solveError_Grid.m
+++ /dev/null
@@ -1,57 +0,0 @@
-\r
-dataIso =[];\r
-dataAniso =[];\r
-maxSize = 500;\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
-%% Anisotrope Verfeinerung\r
-disp 'Anisotrop'\r
-load(Netz)\r
-\r
-ref = 0;\r
-while size(elements,1)<maxSize\r
-  tic\r
-  ref = ref+1; \r
-\r
-  %Netz Isotrop verfeinern (4 Teile)\r
-  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
-  \r
-  %Matrix mit feinem Netz berechnen\r
-  A_fine = mex_build_A(coordinates_fine,elements_fine);\r
-  \r
-  %rechte Seite Aufbauen\r
-  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
-  \r
-  %Lösungsvektor bestimmen\r
-  x_fine = A_fine\b;\r
-  \r
-  %Energienorm^2 bestimmen\r
-  xe_fine = x_fine'*A_fine*x_fine;\r
-\r
-  % Welche Elemente sollen verfeinert werden\r
-  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
-  \r
-  % Wie soll verfeinert werden\r
-  marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
-\r
-  %Daten zur Auswertung zwischenspeichern (testweise?)\r
-  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
-  \r
-  %Sicherheit falls keine Elemente markiert sind\r
-  if(sum(marked~=1)==0)\r
-    disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
-    break;\r
-  end\r
-  \r
-  %Verfeinern\r
-  [coordinates, elements] = refineQuad(coordinates,elements,marked);\r
-  file = ['save' int2str(ref)];\r
-  save ([file], 'coordinates', 'elements')\r
-  toc\r
-end\r
-dataAniso\r
-\r
-clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine\r
-\r

diff --git a/src/test_solveError_Grid1.m b/src/test_solveError_Grid1.m
new file mode 100755 (executable)
index 0000000..fe04315
--- /dev/null
+++ b/src/test_solveError_Grid1.m
@@ -0,0 +1,130 @@
+\r
+mex mex_build_As1.cpp slpRectangle.cpp\r
+\r
+dataIso =[];\r
+dataAniso =[];\r
+maxSize = 500;\r
+fileName = './test_save';\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
+\r
+%% Anisotrope Verfeinerung\r
+if(~exist('d1'))\r
+  disp 'Normal'\r
+  load(Netz)\r
+\r
+  ref = 0;\r
+\r
+  while size(elements,1)<maxSize\r
+    tic\r
+    ref = ref+1; \r
+\r
+\r
+\r
+    %Netz Isotrop verfeinern (4 Teile)\r
+    [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
+\r
+    %Matrix mit feinem Netz berechnen\r
+    A_fine = mex_build_A(coordinates_fine,elements_fine);\r
+\r
+    %rechte Seite Aufbauen\r
+    b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
+\r
+    %Lösungsvektor bestimmen\r
+    x_fine = A_fine\b;\r
+\r
+    %Energienorm^2 bestimmen\r
+    xe_fine = x_fine'*A_fine*x_fine;\r
+\r
+    % Welche Elemente sollen verfeinert werden\r
+    ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
+\r
+    % Wie soll verfeinert werden\r
+    marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
+\r
+    %Daten zur Auswertung zwischenspeichern (testweise?)\r
+    dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
+\r
+    %Sicherheit falls keine Elemente markiert sind\r
+    if(sum(marked~=1)==0)\r
+      disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
+      break;\r
+    end\r
+\r
+    %Verfeinern\r
+    %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
+    file = ['save' int2str(ref)];\r
+    load(file);\r
+    toc\r
+  end\r
+  d1 = dataAniso\r
+end\r
+\r
+%% Anisotrope Verfeinerung\r
+disp 'Optimiert'\r
+load(Netz)\r
+\r
+ref = 0;\r
+while size(elements,1)<maxSize\r
+  tic\r
+  ref = ref+1; \r
+  \r
+\r
+\r
+  %Netz Isotrop verfeinern (4 Teile)\r
+  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
+  \r
+  %Matrix mit feinem Netz berechnen\r
+  A_fine = mex_build_As1(coordinates_fine,elements_fine);\r
+  \r
+  %rechte Seite Aufbauen\r
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
+  \r
+  %Lösungsvektor bestimmen\r
+  x_fine = A_fine\b;\r
+  \r
+  %Energienorm^2 bestimmen\r
+  xe_fine = x_fine'*A_fine*x_fine;\r
+\r
+  % Welche Elemente sollen verfeinert werden\r
+  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
+  \r
+  % Wie soll verfeinert werden\r
+  marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
+\r
+  %Daten zur Auswertung zwischenspeichern (testweise?)\r
+  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
+  \r
+  %Sicherheit falls keine Elemente markiert sind\r
+  if(sum(marked~=1)==0)\r
+    disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
+    break;\r
+  end\r
+  \r
+  %Verfeinern\r
+  %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
+  file = ['save' int2str(ref)];\r
+  load(file);\r
+  toc\r
+end\r
+d2 = dataAniso\r
+\r
+clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine\r
+\r
+\r
+%% Alles Zeichnen\r
+figure (1)\r
+loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));\r
+\r
+legend('Normal','optimiert')\r
+xlabel('Elemente')\r
+ylabel('Fehler')\r
+\r
+figure (2)\r
+loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));\r
+\r
+legend('Normal','Optimiert')\r
+xlabel('Elemente')\r
+ylabel('EnergieNorm')\r

diff --git a/src/test_solveError_Grid2.m b/src/test_solveError_Grid2.m
index f90d5af44d71e45fa973e65ecad8b7062462d51e..52aa416dd4f1087bf7c325daff5f26ff8a57074f 100755 (executable)
--- a/src/test_solveError_Grid2.m
+++ b/src/test_solveError_Grid2.m
@@ -1,4 +1,6 @@
 \r
+mex mex_build_As2.cpp slpRectangle.cpp\r
+\r
 dataIso =[];\r
 dataAniso =[];\r
 maxSize = 500;\r
@@ -8,61 +10,15 @@ Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
 % Netz = 'exmpl_2DQuad';\r
 \r
 \r
-% %% Anisotrope Verfeinerung\r
-% disp 'Anisotrop1'\r
-% load(Netz)\r
-% \r
-% ref = 0;\r
-% while size(elements,1)<maxSize\r
-%   tic\r
-%   ref = ref+1; \r
-%   \r
-% \r
-% \r
-%   %Netz Isotrop verfeinern (4 Teile)\r
-%   [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
-%   \r
-%   %Matrix mit feinem Netz berechnen\r
-%   A_fine = mex_build_A(coordinates_fine,elements_fine);\r
-%   \r
-%   %rechte Seite Aufbauen\r
-%   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
-%   \r
-%   %Lösungsvektor bestimmen\r
-%   x_fine = A_fine\b;\r
-%   \r
-%   %Energienorm^2 bestimmen\r
-%   xe_fine = x_fine'*A_fine*x_fine;\r
-% \r
-%   % Welche Elemente sollen verfeinert werden\r
-%   ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
-%   \r
-%   % Wie soll verfeinert werden\r
-%   marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
-% \r
-%   %Daten zur Auswertung zwischenspeichern (testweise?)\r
-%   dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
-%   \r
-%   %Sicherheit falls keine Elemente markiert sind\r
-%   if(sum(marked~=1)==0)\r
-%     disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
-%     break;\r
-%   end\r
-%   \r
-%   %Verfeinern\r
-%   %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
-%   file = ['save' int2str(ref)];\r
-%   load(file);\r
-%   toc\r
-% end\r
-% d1 = dataAniso\r
+%% Anisotrope Verfeinerung\r
 \r
+if(~exist('d1'))\r
+  \r
+  disp 'Normal'\r
+  load(Netz)\r
 \r
-%% Anisotrope Verfeinerung\r
-disp 'Anisotrop2'\r
-load(Netz)\r
+  ref = 0;\r
 \r
-ref = 0;\r
 while size(elements,1)<maxSize\r
   tic\r
   ref = ref+1; \r
@@ -73,7 +29,7 @@ while size(elements,1)<maxSize
   [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
   \r
   %Matrix mit feinem Netz berechnen\r
-  A_fine = mex_build_As1(coordinates_fine,elements_fine);\r
+  A_fine = mex_build_A(coordinates_fine,elements_fine);\r
   \r
   %rechte Seite Aufbauen\r
   b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
@@ -105,6 +61,56 @@ while size(elements,1)<maxSize
   load(file);\r
   toc\r
 end\r
+d1 = dataAniso\r
+end\r
+\r
+%% Anisotrope Verfeinerung\r
+disp 'Optimiert'\r
+load(Netz)\r
+\r
+ref = 0;\r
+while size(elements,1)<maxSize\r
+  tic\r
+  ref = ref+1; \r
+  \r
+\r
+\r
+  %Netz Isotrop verfeinern (4 Teile)\r
+  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
+  \r
+  %Matrix mit feinem Netz berechnen\r
+  A_fine = mex_build_As2(coordinates_fine,elements_fine);\r
+  \r
+  %rechte Seite Aufbauen\r
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
+  \r
+  %Lösungsvektor bestimmen\r
+  x_fine = A_fine\b;\r
+  \r
+  %Energienorm^2 bestimmen\r
+  xe_fine = x_fine'*A_fine*x_fine;\r
+\r
+  % Welche Elemente sollen verfeinert werden\r
+  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
+  \r
+  % Wie soll verfeinert werden\r
+  marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
+\r
+  %Daten zur Auswertung zwischenspeichern (testweise?)\r
+  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
+  \r
+  %Sicherheit falls keine Elemente markiert sind\r
+%   if(sum(marked~=1)==0)\r
+%     disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
+%     break;\r
+%   end\r
+  \r
+  %Verfeinern\r
+  %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
+  file = ['save' int2str(ref)];\r
+  load(file);\r
+  toc\r
+end\r
 d2 = dataAniso\r
 \r
 clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine\r
@@ -114,13 +120,13 @@ clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine
 figure (1)\r
 loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));\r
 \r
-legend('Isotrop','Anisotrop')\r
+legend('Normal','optimiert')\r
 xlabel('Elemente')\r
 ylabel('Fehler')\r
 \r
 figure (2)\r
 loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));\r
 \r
-legend('Isotrop','Anisotrop')\r
+legend('Normal','Optimiert')\r
 xlabel('Elemente')\r
 ylabel('EnergieNorm')\r