From: Peter Schaefer Date: Wed, 20 Mar 2013 19:46:16 +0000 (+0100) Subject: [src] Quad und Stammfunktion mit Templates durchgezogen (wirkt schneller) X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=ef1c153279bf44e257a616dfa55d21ccc18a6c5b;p=bacc.git [src] Quad und Stammfunktion mit Templates durchgezogen (wirkt schneller) --- diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp index 20c206c..cd8df8f 100644 --- a/src/slpRectangle.cpp +++ b/src/slpRectangle.cpp @@ -145,8 +145,20 @@ double inline dist_s(double b, double t, double d1) { return sqrt(dist_s2(b, t, d1)); } -double inline f_A(double x1, double x2, double y1, double y2, double l) { - return 1. / sqrt((x1 - y1) * (x1 - y1) + (x2 - y2) * (x2 - y2) + l * l); +double inline f_par(double x1, double x2, double y1, double y2, double d1, + double d2, double d3) { + return 1. + / sqrt( + (x1 - y1 - d1) * (x1 - y1 - d1) + + (x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3); +} + +double inline f_ort(double x1, double x2, double y2, double y3, double d1, + double d2, double d3) { + return 1. + / sqrt( + (x1 - d1) * (x1 - d1) + (x2 - y2 - d2) * (x2 - y2 - d2) + + (y3 + d3) * (y3 + d3)); } double inline g_QY(double p, double y, double x, double l) { @@ -154,11 +166,11 @@ double inline g_QY(double p, double y, double x, double l) { double sol = 0; for (int i = 0; i < GAUSS_SIZE[quad]; ++i) { -sol += pow( - (x - GAUSS_NODES[quad][i].n * y) - * (x - GAUSS_NODES[quad][i].n * y) + l * l, p) - * GAUSS_NODES[quad][i].w * y; -} + sol += pow( + (x - GAUSS_NODES[quad][i].n * y) + * (x - GAUSS_NODES[quad][i].n * y) + l * l, p) + * GAUSS_NODES[quad][i].w * y; + } return sol; } @@ -201,15 +213,15 @@ double inline G_QY1Y2(double p, double y1, double y2, double x1, double x2, double sol = 0; for (int i = 0; i < GAUSS_SIZE[quad]; ++i) { for (int j = 0; j < GAUSS_SIZE[quad]; ++j) { -sol += pow( - (x1 - y1 * GAUSS_NODES[quad][i].n) - * (x1 - y1 * GAUSS_NODES[quad][i].n) - + (x2 - y2 * GAUSS_NODES[quad][j].n) - * (x2 - y2 * GAUSS_NODES[quad][j].n) - + l * l, p) * y1 * GAUSS_NODES[quad][i].w * y2 - * GAUSS_NODES[quad][j].w; + sol += pow( + (x1 - y1 * GAUSS_NODES[quad][i].n) + * (x1 - y1 * GAUSS_NODES[quad][i].n) + + (x2 - y2 * GAUSS_NODES[quad][j].n) + * (x2 - y2 * GAUSS_NODES[quad][j].n) + + l * l, p) * y1 * GAUSS_NODES[quad][i].w * y2 + * GAUSS_NODES[quad][j].w; + } } -} return sol; } @@ -225,7 +237,8 @@ double inline G_AY2X2(double y1, double y2, double x1, double x2, double d3) { //sol -= (x2 - y2) * g_AY(-0.5,x2,y2, sqrt((x1 - y1) * (x1 - y1) + d3 * d3)); sol += (x2 - y2) * log(sqrt(hlp) - (x2 - y2)); - if (sol != sol/* || fabs(sol) == std::numeric_limits::infinity()*/) { + if (sol + != sol/* || fabs(sol) == std::numeric_limits::infinity()*/) { // cout << "G_AY2X2: " << sol << " isn't a Number. (" << y1 << "," << y2 // << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl; // cout << (x2 - y2) << "__" << sqrt(hlp) << endl; @@ -279,45 +292,6 @@ double inline Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2, return sol; } -template -double Q2(double a, double b) { - double sol = 0; - - for (int i = 0; i < GAUSS_SIZE[quad]; ++i) { - for (int j = 0; j < GAUSS_SIZE[quad]; ++j) { - sol += f(a * GAUSS_NODES[quad][i].n, b * GAUSS_NODES[quad][j].n) * a - * b * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - } - } - - return sol; -} - -template -double Q4(double a, double b, double c, double d) { - //Fall 0 - - double sol = 0; - int i, j, k, l; - - for (i = 0; i < GAUSS_SIZE[quad]; ++i) { - for (j = 0; j < GAUSS_SIZE[quad]; ++j) { - for (k = 0; k < GAUSS_SIZE[quad]; ++k) { - for (l = 0; l < GAUSS_SIZE[quad]; ++l) { - sol += f(a * GAUSS_NODES[quad][i].n, - b * GAUSS_NODES[quad][j].n, - c * GAUSS_NODES[quad][k].n, - d * GAUSS_NODES[quad][l].n) * a * b * c * d - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w - * GAUSS_NODES[quad][k].w * GAUSS_NODES[quad][l].w; - } - } - } - } - - return sol; -} - double inline F_par(double x1, double x2, double y1, double y2, double d1, double d2, double d3) { @@ -374,160 +348,53 @@ double inline F_ort(double x1, double x2, double y2, double y3, double d1, return sol / 2.; } -/*double applyInt4( - double(*f)(double, double, double, double, double, double, double), - double a, double b, double c, double d, double r, double t, double u, - double v, double d1, double d2, double d3) { - - return f(b, d, t, v, d1, d2, d3) - f(b, d, t, u, d1, d2, d3) - f(b, d, r, - v, d1, d2, d3) + f(b, d, r, u, d1, d2, d3) - f(b, c, t, v, d1, d2, - d3) + f(b, c, t, u, d1, d2, d3) + f(b, c, r, v, d1, d2, d3) - f(b, - c, r, u, d1, d2, d3) - f(a, d, t, v, d1, d2, d3) + f(a, d, t, u, - d1, d2, d3) + f(a, d, r, v, d1, d2, d3) - f(a, d, r, u, d1, d2, d3) - + f(a, c, t, v, d1, d2, d3) - f(a, c, t, u, d1, d2, d3) - f(a, c, - r, v, d1, d2, d3) + f(a, c, r, u, d1, d2, d3); - - }*/ - -template -double inline apply0Int4( - double b, double d, double t, double v, double d1, double d2, - double d3) { - - return f(b, d, t, v, d1, d2, d3) - f(b, d, t, 0, d1, d2, d3) - - f(b, d, 0, v, d1, d2, d3) + f(b, d, 0, 0, d1, d2, d3) - - f(b, 0, t, v, d1, d2, d3) + f(b, 0, t, 0, d1, d2, d3) - + f(b, 0, 0, v, d1, d2, d3) - f(b, 0, 0, 0, d1, d2, d3) - - f(0, d, t, v, d1, d2, d3) + f(0, d, t, 0, d1, d2, d3) - + f(0, d, 0, v, d1, d2, d3) - f(0, d, 0, 0, d1, d2, d3) - + f(0, 0, t, v, d1, d2, d3) - f(0, 0, t, 0, d1, d2, d3) - - f(0, 0, 0, v, d1, d2, d3) + f(0, 0, 0, 0, d1, d2, d3); -} -template -double inline apply0Int2( - double b, double d, double t, double v, double d1, double d2, +double inline FY1Y2_par(double x1, double x2, double y1, double y2, double d1, double d2, double d3) { - - return f(b, d, t, v, d1, d2, d3) - f(b, 0, t, v, d1, d2, d3) - - f(0, d, t, v, d1, d2, d3) + f(0, 0, t, v, d1, d2, d3); - + return G_AY1Y2(-0.5, y1 + d1, y2 + d2, x1, x2, d3); } -// Berechnet das eigentliche Integral fuer parallele Elemente voll Analytisch -double calcParIntA(double b, double d, double t, double v, double d1, double d2, +double inline FY2X2_par(double y1, double x1, double y2, double x2, double d1, double d2, double d3) { - return apply0Int4(b, d, t, v, d1, d2, d3); - + return G_AY2X2(y1 + d1, y2 + d2, x1, x2, d3); } -double calcParIntQX1X2(double b, double d, double t, double v, double d1, - double d2, double d3) { - //Fall 2 - double sol = 0; - for (int i = 0; i < GAUSS_SIZE[quad]; ++i) { - for (int j = 0; j < GAUSS_SIZE[quad]; ++j) { -sol += G_AY1Y2(-0.5, t + d1, v + d2, b * GAUSS_NODES[quad][i].n, - d * GAUSS_NODES[quad][j].n, d3) * b * d - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - sol -= G_AY1Y2(-0.5, d1, v + d2, b * GAUSS_NODES[quad][i].n, - d * GAUSS_NODES[quad][j].n, d3) * b * d - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - sol -= G_AY1Y2(-0.5, t + d1, d2, b * GAUSS_NODES[quad][i].n, - d * GAUSS_NODES[quad][j].n, d3) * b * d - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - sol += G_AY1Y2(-0.5, d1, d2, b * GAUSS_NODES[quad][i].n, - d * GAUSS_NODES[quad][j].n, d3) * b * d - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - } -} - - return sol; -// return Q2<>(b,d); - -} -double calcParIntQY1X1(double b, double d, double t, double v, double d1, - double d2, double d3) { - //Fall 3 +/* + * Quadratur über zwei Integrale und Grenzen in Stammfunktion über zwei Integrale + * b,d - Grenzen für Quadratur + * t,v - Grenzen in Stammfunktion + */ +template +double intQ2A2(double b, double d, double t, double v, double d1, double d2, + double d3) { double sol = 0; for (int i = 0; i < GAUSS_SIZE[quad]; ++i) { for (int j = 0; j < GAUSS_SIZE[quad]; ++j) { - - sol += G_AY2X2(t * GAUSS_NODES[quad][j].n + d1, v + d2, - b * GAUSS_NODES[quad][i].n, d, d3) * t * b - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - sol -= G_AY2X2(t * GAUSS_NODES[quad][j].n + d1, v + d2, - b * GAUSS_NODES[quad][i].n, 0, d3) * t * b - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - sol -= G_AY2X2(t * GAUSS_NODES[quad][j].n + d1, d2, - b * GAUSS_NODES[quad][i].n, d, d3) * t * b - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - sol += G_AY2X2(t * GAUSS_NODES[quad][j].n + d1, d2, - b * GAUSS_NODES[quad][i].n, 0, d3) * t * b - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - - if (sol!= sol /*|| fabs(sol) == numeric_limits::infinity()*/) { - //mexPrintf("(%f,%f),(%f,%f),(%f,%f,%f)",b,d,t,v,d1,d2,d3); - /* cout << "->(" << i << "," << j << ")" << endl; - cout << "-|(" - << G_AY2X2(t * GAUSS_NODES[quad][j].n + d1, v + d2, - b * GAUSS_NODES[quad][i].n, d, d3) * t * b - * GAUSS_NODES[quad][i].w - * GAUSS_NODES[quad][j].w << "," - << G_AY2X2(t * GAUSS_NODES[quad][j].n + d1, v + d2, - b * GAUSS_NODES[quad][i].n, 0, d3) * t * b - * GAUSS_NODES[quad][i].w - * GAUSS_NODES[quad][j].w << "," - << G_AY2X2(t * GAUSS_NODES[quad][j].n + d1, d2, - b * GAUSS_NODES[quad][i].n, d, d3) * t * b - * GAUSS_NODES[quad][i].w - * GAUSS_NODES[quad][j].w << "," - << G_AY2X2(t * GAUSS_NODES[quad][j].n + d1, d2, - b * GAUSS_NODES[quad][i].n, 0, d3) * t * b - * GAUSS_NODES[quad][i].w - * GAUSS_NODES[quad][j].w << ")" << endl; - cout << "||" - << t * b * GAUSS_NODES[quad][i].w - * GAUSS_NODES[quad][j].w - << "||-----------------" << endl; - - cout << "? (" - << G_AY2X2(t * GAUSS_NODES[quad][j].n + d1, d2, - b * GAUSS_NODES[quad][i].n, 0, d3) << "," - << t * b * GAUSS_NODES[quad][i].w - * GAUSS_NODES[quad][j].w << ")" << endl; - - cout << endl;*/ - //mexPrintf("%f\n",sol); - } - return sol; + sol += (F(b * GAUSS_NODES[quad][i].n, d * GAUSS_NODES[quad][j].n, t, + v, d1, d2, d3) + - F(b * GAUSS_NODES[quad][i].n, d * GAUSS_NODES[quad][j].n, + 0, v, d1, d2, d3) + - F(b * GAUSS_NODES[quad][i].n, d * GAUSS_NODES[quad][j].n, + t, 0, d1, d2, d3) + + F(b * GAUSS_NODES[quad][i].n, d * GAUSS_NODES[quad][j].n, + 0, 0, d1, d2, d3)) * GAUSS_NODES[quad][i].w + * GAUSS_NODES[quad][j].w * b * d; } } - //if (sol != sol /*|| sol == numeric_limits::infinity()*/) - // cout << "calcParIntQY1X1: " << sol << " isn't a Number. (" << b << "," - // << d << ")(" << t << "," << v << ")(" << d1 << "," << d2 << "," - // << d3 << ")" << endl; - return sol; } -double calcParIntQY1(double b, double d, double t, double v, double d1, - double d2, double d3) { - //Fall 4 - - double sol = 0; - mexPrintf("warning: calcParIntQY1 nicht implementiert!"); - return 0; ///ACHTUNG noch Falsch - -} - -double calcParIntQ(double b, double d, double t, double v, double d1, double d2, +/* + * Quadratur über vier Integrale + * b,d,t,v - Grenzen für Quadratur + */ +template +double intQ4(double b, double d, double t, double v, double d1, double d2, double d3) { - //Fall 0 double sol = 0; int i, j, k, l; @@ -536,62 +403,48 @@ double calcParIntQ(double b, double d, double t, double v, double d1, double d2, for (j = 0; j < GAUSS_SIZE[quad]; ++j) { for (k = 0; k < GAUSS_SIZE[quad]; ++k) { for (l = 0; l < GAUSS_SIZE[quad]; ++l) { -sol += f_A(b * GAUSS_NODES[quad][i].n, - d * GAUSS_NODES[quad][j].n, - (t) * GAUSS_NODES[quad][k].n + d1, - (v) * GAUSS_NODES[quad][l].n + d2, d3) * b - * GAUSS_NODES[quad][i].w * d - * GAUSS_NODES[quad][j].w * t - * GAUSS_NODES[quad][k].w * v - * GAUSS_NODES[quad][l].w; + sol += F(b * GAUSS_NODES[quad][i].n, + d * GAUSS_NODES[quad][j].n, + t * GAUSS_NODES[quad][k].n, + v * GAUSS_NODES[quad][l].n, d1, d2, d3) * b * d * t + * v * GAUSS_NODES[quad][i].w + * GAUSS_NODES[quad][j].w * GAUSS_NODES[quad][k].w + * GAUSS_NODES[quad][l].w; + } } } } -} return sol; } -double calcOrtIntA(double b, double d, double t, double v, double d1, double d2, - double d3) { - return apply0Int4(b, d, t, v, d1, d2, d3); - -} - -double calcOrtIntQX1X2(double b, double d, double t, double v, double d1, +/* + * Grenzen in Stammfunktion über vier Integrale + * b,d,t,v - Grenzen für Stammfunktion + */ +template +double inline intA4(double b, double d, double t, double v, double d1, double d2, double d3) { - //Fall 2 - double sol = 0; - - for (int i = 0; i < GAUSS_SIZE[quad]; ++i) { - for (int j = 0; j < GAUSS_SIZE[quad]; ++j) { -sol += G_AY1Y2(-0.5, t + d2, v + d3, d * GAUSS_NODES[quad][j].n, 0, - b * GAUSS_NODES[quad][i].n - d1) * b * d - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - sol -= G_AY1Y2(-0.5, d2, v + d3, d * GAUSS_NODES[quad][j].n, 0, - b * GAUSS_NODES[quad][i].n - d1) * b * d - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - sol -= G_AY1Y2(-0.5, t + d2, d3, d * GAUSS_NODES[quad][j].n, 0, - b * GAUSS_NODES[quad][i].n - d1) * b * d - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - sol += G_AY1Y2(-0.5, d2, d3, d * GAUSS_NODES[quad][j].n, 0, - b * GAUSS_NODES[quad][i].n - d1) * b * d - * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w; - } -} - return sol; + return F(b, d, t, v, d1, d2, d3) - F(b, d, t, 0, d1, d2, d3) + - F(b, d, 0, v, d1, d2, d3) + F(b, d, 0, 0, d1, d2, d3) + - F(b, 0, t, v, d1, d2, d3) + F(b, 0, t, 0, d1, d2, d3) + + F(b, 0, 0, v, d1, d2, d3) - F(b, 0, 0, 0, d1, d2, d3) + - F(0, d, t, v, d1, d2, d3) + F(0, d, t, 0, d1, d2, d3) + + F(0, d, 0, v, d1, d2, d3) - F(0, d, 0, 0, d1, d2, d3) + + F(0, 0, t, v, d1, d2, d3) - F(0, 0, t, 0, d1, d2, d3) + - F(0, 0, 0, v, d1, d2, d3) + F(0, 0, 0, 0, d1, d2, d3); } double cParO1(double b, double d, double t, double v, double d1, double d2, double d3, double* zeta) { - return calcParIntA(b, d, t, v, d1, d2, d3); + return intA4(b, d, t, v, d1, d2, d3); } double cOrtO1(double b, double d, double t, double v, double d1, double d2, double d3, double* zeta) { - return calcOrtIntA(b, d, t, v, d1, d2, d3); + return intA4(b, d, t, v, d1, d2, d3); } double cParO2(double b, double d, double t, double v, double d1, double d2, @@ -604,10 +457,10 @@ double cParO2(double b, double d, double t, double v, double d1, double d2, //kurze Achse nach vorn // switch_dim(b, d, t, v, d1, d2); - if (zeta[0]*zeta[0] *(t * t + v * v) < dist2(b, d, t, v, d1, d2, d3)) - return calcParIntQ(b, d, t, v, d1, d2, d3); + if (zeta[0] * zeta[0] * (t * t + v * v) < dist2(b, d, t, v, d1, d2, d3)) + return intQ4(b, d, t, v, d1, d2, d3); - return calcParIntA(b, d, t, v, d1, d2, d3); + return intA4(b, d, t, v, d1, d2, d3); } double cParO3(double b, double d, double t, double v, double d1, double d2, @@ -619,16 +472,13 @@ double cParO3(double b, double d, double t, double v, double d1, double d2, //kurze Achse nach vorn // switch_dim(b, d, t, v, d1, d2); - if (zeta[1]*zeta[1] *(b * b + d * d) < dist2(b, d, t, v, d1, d2, d3)) -// if ((max(b, t) <= zeta[2] * min(d, v)) -// && (min(dist_s2(b, t, d1), dist_s2(d, v, d2)) > 0)) -// mexPrintf("Qx1x2"); - return calcParIntQX1X2(b, d, t, v, d1, d2, d3); + if (zeta[1] * zeta[1] * (b * b + d * d) < dist2(b, d, t, v, d1, d2, d3)) + return intQ2A2(b, d, t, v, d1, d2, d3); - if (zeta[0]*zeta[0] *(t * t + v * v) < dist2(b, d, t, v, d1, d2, d3)) - return calcParIntQ(b, d, t, v, d1, d2, d3); + if (zeta[0] * zeta[0] * (t * t + v * v) < dist2(b, d, t, v, d1, d2, d3)) + return intQ4(b, d, t, v, d1, d2, d3); - return calcParIntA(b, d, t, v, d1, d2, d3); + return intA4(b, d, t, v, d1, d2, d3); } double cParO4(double b, double d, double t, double v, double d1, double d2, @@ -643,10 +493,10 @@ double cParO4(double b, double d, double t, double v, double d1, double d2, switch_dim(b, d, t, v, d1, d2); if ((t * t + v * v) < zeta[0] * dist2(b, d, t, v, d1, d2, d3)) - return calcParIntQ(b, d, t, v, d1, d2, d3); + return intQ4(b, d, t, v, d1, d2, d3); - if ((b * b + d * d) < zeta[0] * dist2(b, d, t, v, d1, d2, d3)) - return calcParIntQX1X2(b, d, t, v, d1, d2, d3); +// if ((b * b + d * d) < zeta[0] * dist2(b, d, t, v, d1, d2, d3)) +// return calcParIntQX1X2(b, d, t, v, d1, d2, d3); - return calcParIntA(b, d, t, v, d1, d2, d3); + return intA4(b, d, t, v, d1, d2, d3); } diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp index 45f22f3..2559533 100644 --- a/src/slpRectangle.hpp +++ b/src/slpRectangle.hpp @@ -42,15 +42,15 @@ //double, double, double, double, double, double, double); // sol = calcParInt.(b,d,t,v,d1,d2,d3); -double calcParIntA(double, double, double, double, double, double, double); -double calcParIntQX1X2(double, double, double, double, double, double, double); -double calcParIntQY1X1(double, double, double, double, double, double, double); -double calcParIntQY1(double, double, double, double, double, double, double); -double calcParIntQ(double, double, double, double, double, double, double); +//double calcParIntA(double, double, double, double, double, double, double); +//double calcParIntQX1X2(double, double, double, double, double, double, double); +//double calcParIntQY1X1(double, double, double, double, double, double, double); +//double calcParIntQY1(double, double, double, double, double, double, double); +//double calcParIntQ(double, double, double, double, double, double, double); // sol = calcOrtInt.(b,d,t,v,d1,d2,d3); -double calcOrtIntA(double, double, double, double, double, double, double); -double calcOrtIntQX1X2(double, double, double, double, double, double, double); +//double calcOrtIntA(double, double, double, double, double, double, double); +//double calcOrtIntQX1X2(double, double, double, double, double, double, double); //Voll Analytisch double cParO1(double, double, double, double, double, double, double, double*);