[src] test_v

[src] Beschreibung in mex_build_V
[src] zeta in Type 0 gedreht

diff --git a/src/A_step.m b/src/A_step.m
index 096c33b334e039a3542d7f5240f06ca6aee47d75..0e208ce15a45419fdda6b16f0726785b28188f74 100644 (file)
--- a/src/A_step.m
+++ b/src/A_step.m
@@ -39,6 +39,10 @@ time = zeros(1,3);
   %Flaecheninhalte Berechnen (rhs)
 %   b = sqrt(sum(quadNorm(coo_fine,ele_fine,'w').^2,2));
   b_fine = areaQuad(coo_fine,ele_fine);
+  
+  [b,tnormal,s] = areaQuad(G_C,G_E);
+  hmin = min(s,[],2);
+  hmax = max(s,[],2);
 
   tic
   %data -> ErgebnisMatrix
@@ -46,6 +50,8 @@ time = zeros(1,3);
   %save_* -> ZwischenSpeicherung
  save_A = {};
  save_x = {};
+ save_A_fine = {};
+ save_x_fine = {};
   
   %Alle MatrixBrechenungsArten mit dem selben Netz berechnen
   for i = 1:length(type)
@@ -63,10 +69,6 @@ time = zeros(1,3);
       
     %Lösung Berechnen
     x_fine = A_fine\b_fine;
-
-    [b,tnormal,s] = areaQuad(G_C,G_E);
-    hmin = min(s,[],2);
-    hmax = max(s,[],2);
     
     % \tilde \mu ( \Pi h -h + L_2 )
     tmu = hmin.* b .* sum((x_fine(f2s)'-repmat(sum(x_fine(f2s)',1)/4,4,1)).^2)' /4;
@@ -94,8 +96,12 @@ time = zeros(1,3);
     
     eta = abs(xe_fine-xe);
     
-   save_A{i} = A_fine;
-   save_x{i} = x_fine;
+   save_A_fine{i} = A_fine;
+   save_x_fine{i} = x_fine;
+   
+   save_A{i} = A;
+   save_x{i} = x;
+   
     
     
     
@@ -146,7 +152,8 @@ time = zeros(1,3);
    end
   typeN = int2str(type);
  save (['meshSave/fine_' out typeN(typeN~=' ') '_' int2str(size(G_T,1))],...
-     'coo_fine', 'ele_fine','neigh_fine','f2s','time','data','save_A','save_x')
+     'coo_fine', 'ele_fine','neigh_fine','f2s','data',...
+     'save_A','save_x','save_A_fine','save_x_fine','b','b_fine')
 
 %   clear 'coo_fine' 'ele_fine' 'neigh_fine' 'f2s' 
 %   plotShape(G_C,G_E,'');

diff --git a/src/mex_build_V.cpp b/src/mex_build_V.cpp
index 6f11e60791f1c525b59c2708301e8ba276464a55..27fcada6e58c8280e3da97b2fadd06778c093f8a 100644 (file)
--- a/src/mex_build_V.cpp
+++ b/src/mex_build_V.cpp
@@ -1,8 +1,26 @@
 /*\r
- * Art der Berechnung durch Parameter bestimmt...\r
- * Analytisch/Semianalytisch/...\r
+ * MEX FUNKTION mex_build_V(coo,ele,zeta,type)\r
  *\r
+ * coo  : Koordinaten Matrix\r
+ * ele  : Element Matrix\r
+ * zeta : Steuerungsvariable der Zulaessigeitsbedingungen\r
+ * type : Bestimmt die Art der Berechnung\r
  *\r
+ * TYPES----\r
+ * 1 : Voll Analytisch\r
+ *             zeta wird ignoriert\r
+ *\r
+ * 0 : Semianalytisch im Fernfeld Quadratur ueber beide Elemente\r
+ *             Quadratur bei : ((t * t + v * v) < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) : wobei v,t jeweils die laengste Seite\r
+ *                              in x und in y Richtung ist und d(1,2,3) der Abstandsvektor zwischen den Elementen\r
+ *\r
+ * 2 : Semianalytisch im Fernfeld Quadratur ueber ein Element\r
+ *             Quadratur bei : ((b * b + d * d) < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) : wobei b,d jeweils die kuerzeste Seite\r
+ *                              in x und in y Richtung ist und d(1,2,3) der Abstandsvektor zwischen den Elementen\r
+ *\r
+ *\r
+ *\r
+ * Peter Schaefer\r
  */\r
 //#include <iostream>\r
 #include <cmath>\r
@@ -10,17 +28,17 @@
 #include "mex.h"\r
 /*/#include <stdlib.h>\r
 \r
-//#include "gauss.hpp"\r
-*/\r
+ //#include "gauss.hpp"\r
+ */\r
 //#define M_EPS 10e-8\r
 /*\r
-//#include "tbb/parallel_for.h"\r
-*/\r
+ //#include "tbb/parallel_for.h"\r
+ */\r
 #include "slpRectangle.hpp"\r
 /*\r
-//using namespace std;\r
-//using namespace slpR;\r
-*/\r
+ //using namespace std;\r
+ //using namespace slpR;\r
+ */\r
 int dimOfVec(double* vec) {\r
        if (vec[2] != 0)\r
                return 2;\r
@@ -31,8 +49,8 @@ int dimOfVec(double* vec) {
        else {\r
                mexErrMsgTxt("length of Site is zero");\r
                /*\r
-//             cerr << "dimOfVec : (" << vec[0] << " " << vec[1] << " " << vec[2]\r
-//                             << ") alle Werte 0 \n";*/\r
+                //             cerr << "dimOfVec : (" << vec[0] << " " << vec[1] << " " << vec[2]\r
+                //                             << ") alle Werte 0 \n";*/\r
                return -1;\r
        }\r
 \r
@@ -44,10 +62,10 @@ inline int dimOfThird(int a, int b) {
 }\r
 \r
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
-/*\r
-//     initeQuad();\r
-//     cout << Quad[1].nod[0];\r
-*/\r
+       /*\r
+        //     initeQuad();\r
+        //     cout << Quad[1].nod[0];\r
+        */\r
        int i, j, k; //Schleifenindizes\r
        double tmp; //Zwischenspeicherung der Einzelnen Werte\r
        int count;\r
@@ -101,9 +119,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
        double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
 \r
        // Welche Funktion soll verwendet werden\r
-       double(*ctypeP)(double, double, double, double, double, double, double,\r
+       double (*ctypeP)(double, double, double, double, double, double, double,\r
                        double*);\r
-       double(*ctypeO)(double, double, double, double, double, double, double,\r
+       double (*ctypeO)(double, double, double, double, double, double, double,\r
                        double*);\r
 \r
        //Art der Berechnung bestimmen\r
@@ -261,30 +279,30 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                        A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
                        if (k != j)\r
                                A[(j * em) + k] = 1. / (4 * M_PI) * tmp;\r
-/*                     if (count++ > ((em * (em + 1)) / 2) / 10) {\r
-                               count = 0;\r
-//                             mexPrintf("#");\r
-//                             cout << "#";\r
-//                             cout.flush();\r
-                       }*/\r
+                       /*                      if (count++ > ((em * (em + 1)) / 2) / 10) {\r
+                        count = 0;\r
+                        //                             mexPrintf("#");\r
+                        //                             cout << "#";\r
+                        //                             cout.flush();\r
+                        }*/\r
 \r
                }\r
 \r
        }\r
 //     cout << endl;\r
        //Rueckgabe (eventuell zurueck schreiben)\r
-/*     mxFree(x0);\r
-       mxFree(x1);\r
-       mxFree(x3);\r
-       mxFree(xa);\r
-       mxFree(xb);\r
-       mxFree(y0);\r
-       mxFree(y1);\r
-       mxFree(y3);\r
-       mxFree(ya);\r
-       mxFree(yb);\r
-       mxFree(d);\r
-       mxFree(dt);\r
-*/\r
+       /*      mxFree(x0);\r
+        mxFree(x1);\r
+        mxFree(x3);\r
+        mxFree(xa);\r
+        mxFree(xb);\r
+        mxFree(y0);\r
+        mxFree(y1);\r
+        mxFree(y3);\r
+        mxFree(ya);\r
+        mxFree(yb);\r
+        mxFree(d);\r
+        mxFree(dt);\r
+        */\r
        return;\r
 }\r

diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp
index 4937d0476489d0168d72310984cc985087c9e639..6bb6781cafe47246de5ba0fc2215ce6cfc6636e8 100644 (file)
--- a/src/slpRectangle.cpp
+++ b/src/slpRectangle.cpp
@@ -10,7 +10,6 @@
 #include "gauss.hpp"\r
 \r
 #define quad 4 // Anzahl der Quadratur Auswertungstellen 2^quad\r
-\r
 //using namespace std;\r
 \r
 int sign(double);\r
@@ -116,10 +115,10 @@ double inline G_AY2X2(double y1, double y2, double x1, double x2, double d3) {
                sol += (x2 - y2) * log(sqrt(hlp) - (x2 - y2));\r
 \r
        if (sol != sol /*|| fabs(sol) == numeric_limits<double>::infinity()*/) {\r
-       //      cout << "G_AY2X2: " << sol << " isn't a Number. (" << y1 << "," << y2\r
-       //                      << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl;\r
-       //      cout << (x2 - y2) << "__" << sqrt(hlp) << endl;\r
-       //      cout << (x1 - y1) * (x1 - y1) << " " << hlp << endl;\r
+               //      cout << "G_AY2X2: " << sol << " isn't a Number. (" << y1 << "," << y2\r
+               //                      << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl;\r
+               //      cout << (x2 - y2) << "__" << sqrt(hlp) << endl;\r
+               //      cout << (x1 - y1) * (x1 - y1) << " " << hlp << endl;\r
        }\r
        return sol;\r
 }\r
@@ -154,8 +153,8 @@ double inline G_AY1Y2(double p, double y1, double y2, double x1, double x2,
                sol /= 2 * p + 2;\r
        } else {\r
                sol = NAN;\r
-       //      cout << "warning in G_AY1Y2: no case for p=" << p\r
-       //                      << " defined. Possible: [-1.5,-0.5,0.5,...]" << endl;\r
+               //      cout << "warning in G_AY1Y2: no case for p=" << p\r
+               //                      << " defined. Possible: [-1.5,-0.5,0.5,...]" << endl;\r
        }\r
 \r
        return sol;\r
@@ -205,10 +204,7 @@ double inline F_ort(double x1, double x2, double y2, double y3, double d1,
 \r
        if ((x1 - d1) != 0)\r
                sol += (x1 - d1)\r
-                               * g_AY(\r
-                                               0.5,\r
-                                               y3,\r
-                                               -d3,\r
+                               * g_AY(0.5, y3, -d3,\r
                                                sqrt(\r
                                                                (x1 - d1) * (x1 - d1)\r
                                                                                + (x2 - y2 - d2) * (x2 - y2 - d2)));\r
@@ -219,10 +215,7 @@ double inline F_ort(double x1, double x2, double y2, double y3, double d1,
 \r
        if ((y3 + d3) != 0)\r
                sol += (y3 + d3)\r
-                               * g_AY(\r
-                                               0.5,\r
-                                               x1,\r
-                                               d1,\r
+                               * g_AY(0.5, x1, d1,\r
                                                sqrt(\r
                                                                (x2 - y2 - d2) * (x2 - y2 - d2)\r
                                                                                + (y3 + d3) * (y3 + d3)));\r
@@ -246,7 +239,7 @@ double inline F_ort(double x1, double x2, double y2, double y3, double d1,
  }*/\r
 \r
 double inline apply0Int4(\r
-               double(*f)(double, double, double, double, double, double, double),\r
+               double (*f)(double, double, double, double, double, double, double),\r
                double b, double d, double t, double v, double d1, double d2,\r
                double d3) {\r
 \r
@@ -262,7 +255,7 @@ double inline apply0Int4(
 }\r
 \r
 double inline apply0Int2(\r
-               double(*f)(double, double, double, double, double, double, double),\r
+               double (*f)(double, double, double, double, double, double, double),\r
                double b, double d, double t, double v, double d1, double d2,\r
                double d3) {\r
 \r
@@ -272,8 +265,8 @@ double inline apply0Int2(
 }\r
 \r
 // Berechnet das eigentliche Integral fuer parallele Elemente voll Analytisch\r
-double calcParIntA(double b, double d, double t, double v, double d1,\r
-               double d2, double d3) {\r
+double calcParIntA(double b, double d, double t, double v, double d1, double d2,\r
+               double d3) {\r
        return apply0Int4(F_par, b, d, t, v, d1, d2, d3);\r
 \r
 }\r
@@ -325,37 +318,38 @@ double calcParIntQY1X1(double b, double d, double t, double v, double d1,
                                        b * gauss_nodes[quad][i].n, 0, d3) * t * b\r
                                        * gauss_nodes[quad][i].w * gauss_nodes[quad][j].w;\r
 \r
-                       if (sol != sol /*|| fabs(sol) == numeric_limits<double>::infinity()*/) {\r
-/*                             cout << "->(" << i << "," << j << ")" << endl;\r
-                               cout << "-|("\r
-                                               << G_AY2X2(t * gauss_nodes[quad][j].n + d1, v + d2,\r
-                                                               b * gauss_nodes[quad][i].n, d, d3) * t * b\r
-                                                               * gauss_nodes[quad][i].w\r
-                                                               * gauss_nodes[quad][j].w << ","\r
-                                               << G_AY2X2(t * gauss_nodes[quad][j].n + d1, v + d2,\r
-                                                               b * gauss_nodes[quad][i].n, 0, d3) * t * b\r
-                                                               * gauss_nodes[quad][i].w\r
-                                                               * gauss_nodes[quad][j].w << ","\r
-                                               << G_AY2X2(t * gauss_nodes[quad][j].n + d1, d2,\r
-                                                               b * gauss_nodes[quad][i].n, d, d3) * t * b\r
-                                                               * gauss_nodes[quad][i].w\r
-                                                               * gauss_nodes[quad][j].w << ","\r
-                                               << G_AY2X2(t * gauss_nodes[quad][j].n + d1, d2,\r
-                                                               b * gauss_nodes[quad][i].n, 0, d3) * t * b\r
-                                                               * gauss_nodes[quad][i].w\r
-                                                               * gauss_nodes[quad][j].w << ")" << endl;\r
-                               cout << "||"\r
-                                               << t * b * gauss_nodes[quad][i].w\r
-                                                               * gauss_nodes[quad][j].w\r
-                                               << "||-----------------" << endl;\r
-\r
-                               cout << "? ("\r
-                                               << G_AY2X2(t * gauss_nodes[quad][j].n + d1, d2,\r
-                                                               b * gauss_nodes[quad][i].n, 0, d3) << ","\r
-                                               << t * b * gauss_nodes[quad][i].w\r
-                                                               * gauss_nodes[quad][j].w << ")" << endl;\r
-\r
-                               cout << endl;*/\r
+                       if (sol\r
+                                       != sol /*|| fabs(sol) == numeric_limits<double>::infinity()*/) {\r
+                               /*                              cout << "->(" << i << "," << j << ")" << endl;\r
+                                cout << "-|("\r
+                                << G_AY2X2(t * gauss_nodes[quad][j].n + d1, v + d2,\r
+                                b * gauss_nodes[quad][i].n, d, d3) * t * b\r
+                                * gauss_nodes[quad][i].w\r
+                                * gauss_nodes[quad][j].w << ","\r
+                                << G_AY2X2(t * gauss_nodes[quad][j].n + d1, v + d2,\r
+                                b * gauss_nodes[quad][i].n, 0, d3) * t * b\r
+                                * gauss_nodes[quad][i].w\r
+                                * gauss_nodes[quad][j].w << ","\r
+                                << G_AY2X2(t * gauss_nodes[quad][j].n + d1, d2,\r
+                                b * gauss_nodes[quad][i].n, d, d3) * t * b\r
+                                * gauss_nodes[quad][i].w\r
+                                * gauss_nodes[quad][j].w << ","\r
+                                << G_AY2X2(t * gauss_nodes[quad][j].n + d1, d2,\r
+                                b * gauss_nodes[quad][i].n, 0, d3) * t * b\r
+                                * gauss_nodes[quad][i].w\r
+                                * gauss_nodes[quad][j].w << ")" << endl;\r
+                                cout << "||"\r
+                                << t * b * gauss_nodes[quad][i].w\r
+                                * gauss_nodes[quad][j].w\r
+                                << "||-----------------" << endl;\r
+\r
+                                cout << "? ("\r
+                                << G_AY2X2(t * gauss_nodes[quad][j].n + d1, d2,\r
+                                b * gauss_nodes[quad][i].n, 0, d3) << ","\r
+                                << t * b * gauss_nodes[quad][i].w\r
+                                * gauss_nodes[quad][j].w << ")" << endl;\r
+\r
+                                cout << endl;*/\r
                                return sol;\r
                        }\r
 \r
@@ -380,12 +374,12 @@ double calcParIntQY1(double b, double d, double t, double v, double d1,
 \r
 }\r
 \r
-double calcParIntQ(double b, double d, double t, double v, double d1,\r
-               double d2, double d3) {\r
+double calcParIntQ(double b, double d, double t, double v, double d1, double d2,\r
+               double d3) {\r
        //Fall 0\r
 \r
        double sol = 0;\r
-       int i,j,k,l;\r
+       int i, j, k, l;\r
 \r
        for (i = 0; i < gauss_size[quad]; ++i) {\r
                for (j = 0; j < gauss_size[quad]; ++j) {\r
@@ -405,8 +399,8 @@ double calcParIntQ(double b, double d, double t, double v, double d1,
        return sol;\r
 }\r
 \r
-double calcOrtIntA(double b, double d, double t, double v, double d1,\r
-               double d2, double d3) {\r
+double calcOrtIntA(double b, double d, double t, double v, double d1, double d2,\r
+               double d3) {\r
        return apply0Int4(F_ort, b, d, t, v, d1, d2, d3);\r
 \r
 }\r
@@ -540,7 +534,7 @@ double calcOrtIntO2(double b, double d, double t, double v, double d1,
        if ((b * b + d * d) * zeta <= d1 * d1 + d2 * d2 + d3 * d3) {\r
                return calcOrtIntQX1X2(b, d, t, v, d1, d2, d3);\r
        } else {\r
-       //      cout << "2error";\r
+               //      cout << "2error";\r
                return calcOrtIntA(b, d, t, v, d1, d2, d3);\r
        }\r
 }\r
@@ -583,8 +577,8 @@ double calcOrtIntO3(double b, double d, double t, double v, double d1,
        }\r
 }\r
 \r
-double cParO0(double b, double d, double t, double v, double d1,\r
-               double d2, double d3, double* zeta) {\r
+double cParO0(double b, double d, double t, double v, double d1, double d2,\r
+               double d3, double* zeta) {\r
        double tmp = 0;\r
 \r
        // kleine Seite nach vorn\r
@@ -603,7 +597,7 @@ double cParO0(double b, double d, double t, double v, double d1,
                d2 = -d2;\r
        }\r
 \r
-       if ((t * t + v * v) * zeta[0] < d1 * d1 + d2 * d2 + d3 * d3) {\r
+       if ((t * t + v * v) < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) {\r
                return calcParIntQ(b, d, t, v, d1, d2, d3);\r
        } else {\r
                return calcParIntA(b, d, t, v, d1, d2, d3);\r
@@ -612,13 +606,13 @@ double cParO0(double b, double d, double t, double v, double d1,
        return 0;\r
 }\r
 \r
-double cParO1(double b, double d, double t, double v, double d1,\r
-               double d2, double d3, double* zeta) {\r
+double cParO1(double b, double d, double t, double v, double d1, double d2,\r
+               double d3, double* zeta) {\r
        return calcParIntA(b, d, t, v, d1, d2, d3);\r
 }\r
 \r
-double cParO2(double b, double d, double t, double v, double d1,\r
-               double d2, double d3, double* zeta) {\r
+double cParO2(double b, double d, double t, double v, double d1, double d2,\r
+               double d3, double* zeta) {\r
 \r
        double tmp = 0;\r
 \r
@@ -638,7 +632,7 @@ double cParO2(double b, double d, double t, double v, double d1,
                d2 = -d2;\r
        }\r
 \r
-       if ((b * b + d * d)  < zeta[0]* (d1 * d1 + d2 * d2 + d3 * d3)) {\r
+       if ((b * b + d * d) < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) {\r
 //             cout << "E";\r
                return calcParIntQX1X2(b, d, t, v, d1, d2, d3);\r
        } else {\r
@@ -648,8 +642,8 @@ double cParO2(double b, double d, double t, double v, double d1,
        return 0;\r
 }\r
 \r
-double cParO3(double b, double d, double t, double v, double d1,\r
-               double d2, double d3, double* zeta) {\r
+double cParO3(double b, double d, double t, double v, double d1, double d2,\r
+               double d3, double* zeta) {\r
 \r
        double tmp = 0;\r
 \r
@@ -669,7 +663,7 @@ double cParO3(double b, double d, double t, double v, double d1,
                d2 = -d2;\r
        }\r
 \r
-       if ((b * b + d * d)  < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) {\r
+       if ((b * b + d * d) < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) {\r
 \r
                //kleinere Achse nach vorn\r
                if (b * b + t * t > d * d + v * v) {\r
@@ -687,7 +681,7 @@ double cParO3(double b, double d, double t, double v, double d1,
                        d2 = tmp;\r
                }\r
 \r
-               if (min(b, t)  <= zeta[1] * d1) {\r
+               if (min(b, t) <= zeta[1] * d1) {\r
 //                     cout << "A";\r
                        return calcParIntQY1X1(b, d, t, v, d1, d2, d3);\r
                } else {\r

diff --git a/src/test_v.m b/src/test_v.m
new file mode 100644 (file)
index 0000000..4f1cb4a
--- /dev/null
+++ b/src/test_v.m
@@ -0,0 +1,16 @@
+
+%lade Netz
+    load meshSave/test1_1_10
+
+%Berechne V 2 mal
+    %Voll Analytisch
+%     V1 = mex_build_V(coordinates,elements,[0.7],[1]);
+    %SemiAnalytisch
+%     V2 = mex_build_V(coordinates,elements,[0.7],[0]);
+      V2 = build_A2(coordinates,elements);
+
+
+%Vergleiche beide Matritzen
+   Vdif = abs(V1-V2);
+   spy(Vdif)
+   max(max(Vdif))
\ No newline at end of file