return sqrt(dist_s2(b, t, d1));\r
}\r
\r
-double inline f_A(double x1, double x2, double y1, double y2, double l) {\r
- return 1. / sqrt((x1 - y1) * (x1 - y1) + (x2 - y2) * (x2 - y2) + l * l);\r
+double inline f_par(double x1, double x2, double y1, double y2, double d1,\r
+ double d2, double d3) {\r
+ return 1.\r
+ / sqrt(\r
+ (x1 - y1 - d1) * (x1 - y1 - d1)\r
+ + (x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3);\r
+}\r
+\r
+double inline f_ort(double x1, double x2, double y2, double y3, double d1,\r
+ double d2, double d3) {\r
+ return 1.\r
+ / sqrt(\r
+ (x1 - d1) * (x1 - d1) + (x2 - y2 - d2) * (x2 - y2 - d2)\r
+ + (y3 + d3) * (y3 + d3));\r
}\r
\r
double inline g_QY(double p, double y, double x, double l) {\r
double sol = 0;\r
\r
for (int i = 0; i < GAUSS_SIZE[quad]; ++i) {\r
-sol += pow(\r
- (x - GAUSS_NODES[quad][i].n * y)\r
- * (x - GAUSS_NODES[quad][i].n * y) + l * l, p)\r
- * GAUSS_NODES[quad][i].w * y;\r
-}\r
+ sol += pow(\r
+ (x - GAUSS_NODES[quad][i].n * y)\r
+ * (x - GAUSS_NODES[quad][i].n * y) + l * l, p)\r
+ * GAUSS_NODES[quad][i].w * y;\r
+ }\r
\r
return sol;\r
}\r
double sol = 0;\r
for (int i = 0; i < GAUSS_SIZE[quad]; ++i) {\r
for (int j = 0; j < GAUSS_SIZE[quad]; ++j) {\r
-sol += pow(\r
- (x1 - y1 * GAUSS_NODES[quad][i].n)\r
- * (x1 - y1 * GAUSS_NODES[quad][i].n)\r
- + (x2 - y2 * GAUSS_NODES[quad][j].n)\r
- * (x2 - y2 * GAUSS_NODES[quad][j].n)\r
- + l * l, p) * y1 * GAUSS_NODES[quad][i].w * y2\r
- * GAUSS_NODES[quad][j].w;\r
+ sol += pow(\r
+ (x1 - y1 * GAUSS_NODES[quad][i].n)\r
+ * (x1 - y1 * GAUSS_NODES[quad][i].n)\r
+ + (x2 - y2 * GAUSS_NODES[quad][j].n)\r
+ * (x2 - y2 * GAUSS_NODES[quad][j].n)\r
+ + l * l, p) * y1 * GAUSS_NODES[quad][i].w * y2\r
+ * GAUSS_NODES[quad][j].w;\r
+ }\r
}\r
-}\r
return sol;\r
}\r
\r
//sol -= (x2 - y2) * g_AY(-0.5,x2,y2, sqrt((x1 - y1) * (x1 - y1) + d3 * d3));\r
sol += (x2 - y2) * log(sqrt(hlp) - (x2 - y2));\r
\r
- if (sol != sol/* || fabs(sol) == std::numeric_limits<double>::infinity()*/) {\r
+ if (sol\r
+ != sol/* || fabs(sol) == std::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
return sol;\r
}\r
\r
-template<double (f)(double, double)>\r
-double Q2(double a, double b) {\r
- double sol = 0;\r
-\r
- for (int i = 0; i < GAUSS_SIZE[quad]; ++i) {\r
- for (int j = 0; j < GAUSS_SIZE[quad]; ++j) {\r
- sol += f(a * GAUSS_NODES[quad][i].n, b * GAUSS_NODES[quad][j].n) * a\r
- * b * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w;\r
- }\r
- }\r
-\r
- return sol;\r
-}\r
-\r
-template<double (f)(double, double, double, double)>\r
-double Q4(double a, double b, double c, double d) {\r
- //Fall 0\r
-\r
- double sol = 0;\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
- for (k = 0; k < GAUSS_SIZE[quad]; ++k) {\r
- for (l = 0; l < GAUSS_SIZE[quad]; ++l) {\r
- sol += f(a * GAUSS_NODES[quad][i].n,\r
- b * GAUSS_NODES[quad][j].n,\r
- c * GAUSS_NODES[quad][k].n,\r
- d * GAUSS_NODES[quad][l].n) * a * b * c * d\r
- * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w\r
- * GAUSS_NODES[quad][k].w * GAUSS_NODES[quad][l].w;\r
- }\r
- }\r
- }\r
- }\r
-\r
- return sol;\r
-}\r
-\r
double inline F_par(double x1, double x2, double y1, double y2, double d1,\r
double d2, double d3) {\r
\r
return sol / 2.;\r
}\r
\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
-\r
- return f(b, d, t, v, d1, d2, d3) - f(b, d, t, u, d1, d2, d3) - f(b, d, r,\r
- v, d1, d2, d3) + f(b, d, r, u, d1, d2, d3) - f(b, c, t, v, d1, d2,\r
- d3) + f(b, c, t, u, d1, d2, d3) + f(b, c, r, v, d1, d2, d3) - f(b,\r
- c, r, u, d1, d2, d3) - f(a, d, t, v, d1, d2, d3) + f(a, d, t, u,\r
- d1, d2, d3) + f(a, d, r, v, d1, d2, d3) - f(a, d, r, u, d1, d2, d3)\r
- + f(a, c, t, v, d1, d2, d3) - f(a, c, t, u, d1, d2, d3) - f(a, c,\r
- r, v, d1, d2, d3) + f(a, c, r, u, d1, d2, d3);\r
-\r
- }*/\r
-\r
-template <double (f)(double, double, double, double, double, double, double)>\r
-double inline apply0Int4(\r
- double b, double d, double t, double v, double d1, double d2,\r
- double d3) {\r
-\r
- return f(b, d, t, v, d1, d2, d3) - f(b, d, t, 0, d1, d2, d3)\r
- - f(b, d, 0, v, d1, d2, d3) + f(b, d, 0, 0, d1, d2, d3)\r
- - f(b, 0, t, v, d1, d2, d3) + f(b, 0, t, 0, d1, d2, d3)\r
- + f(b, 0, 0, v, d1, d2, d3) - f(b, 0, 0, 0, d1, d2, d3)\r
- - f(0, d, t, v, d1, d2, d3) + f(0, d, t, 0, d1, d2, d3)\r
- + f(0, d, 0, v, d1, d2, d3) - f(0, d, 0, 0, d1, d2, d3)\r
- + f(0, 0, t, v, d1, d2, d3) - f(0, 0, t, 0, d1, d2, d3)\r
- - f(0, 0, 0, v, d1, d2, d3) + f(0, 0, 0, 0, d1, d2, d3);\r
\r
-}\r
-template <double (f)(double, double, double, double, double, double, double)>\r
-double inline apply0Int2(\r
- double b, double d, double t, double v, double d1, double d2,\r
+double inline FY1Y2_par(double x1, double x2, double y1, double y2, double d1, double d2,\r
double d3) {\r
-\r
- return f(b, d, t, v, d1, d2, d3) - f(b, 0, t, v, d1, d2, d3)\r
- - f(0, d, t, v, d1, d2, d3) + f(0, 0, t, v, d1, d2, d3);\r
-\r
+ return G_AY1Y2(-0.5, y1 + d1, y2 + d2, x1, x2, d3);\r
}\r
\r
-// Berechnet das eigentliche Integral fuer parallele Elemente voll Analytisch\r
-double calcParIntA(double b, double d, double t, double v, double d1, double d2,\r
+double inline FY2X2_par(double y1, double x1, double y2, double x2, double d1, double d2,\r
double d3) {\r
- return apply0Int4<F_par>(b, d, t, v, d1, d2, d3);\r
-\r
+ return G_AY2X2(y1 + d1, y2 + d2, x1, x2, d3);\r
}\r
\r
-double calcParIntQX1X2(double b, double d, double t, double v, double d1,\r
- double d2, double d3) {\r
- //Fall 2\r
- double sol = 0;\r
\r
- for (int i = 0; i < GAUSS_SIZE[quad]; ++i) {\r
- for (int j = 0; j < GAUSS_SIZE[quad]; ++j) {\r
-sol += G_AY1Y2(-0.5, t + d1, v + d2, b * GAUSS_NODES[quad][i].n,\r
- d * GAUSS_NODES[quad][j].n, d3) * b * d\r
- * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w;\r
- sol -= G_AY1Y2(-0.5, d1, v + d2, b * GAUSS_NODES[quad][i].n,\r
- d * GAUSS_NODES[quad][j].n, d3) * b * d\r
- * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w;\r
- sol -= G_AY1Y2(-0.5, t + d1, d2, b * GAUSS_NODES[quad][i].n,\r
- d * GAUSS_NODES[quad][j].n, d3) * b * d\r
- * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w;\r
- sol += G_AY1Y2(-0.5, d1, d2, b * GAUSS_NODES[quad][i].n,\r
- d * GAUSS_NODES[quad][j].n, d3) * b * d\r
- * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w;\r
- }\r
-}\r
-\r
- return sol;\r
-// return Q2<>(b,d);\r
-\r
-}\r
\r
-double calcParIntQY1X1(double b, double d, double t, double v, double d1,\r
- double d2, double d3) {\r
- //Fall 3\r
+/*\r
+ * Quadratur über zwei Integrale und Grenzen in Stammfunktion über zwei Integrale\r
+ * b,d - Grenzen für Quadratur\r
+ * t,v - Grenzen in Stammfunktion\r
+ */\r
+template<double (F)(double, double, double, double, double, double, double)>\r
+double intQ2A2(double b, double d, double t, double v, double d1, double d2,\r
+ double d3) {\r
double sol = 0;\r
\r
for (int i = 0; i < GAUSS_SIZE[quad]; ++i) {\r
for (int j = 0; j < GAUSS_SIZE[quad]; ++j) {\r
-\r
- sol += 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 * GAUSS_NODES[quad][j].w;\r
- sol -= 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 * GAUSS_NODES[quad][j].w;\r
- sol -= 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 * GAUSS_NODES[quad][j].w;\r
- sol += 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 * GAUSS_NODES[quad][j].w;\r
-\r
- if (sol!= sol /*|| fabs(sol) == numeric_limits<double>::infinity()*/) {\r
- //mexPrintf("(%f,%f),(%f,%f),(%f,%f,%f)",b,d,t,v,d1,d2,d3);\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
- //mexPrintf("%f\n",sol);\r
- }\r
- return sol;\r
+ sol += (F(b * GAUSS_NODES[quad][i].n, d * GAUSS_NODES[quad][j].n, t,\r
+ v, d1, d2, d3)\r
+ - F(b * GAUSS_NODES[quad][i].n, d * GAUSS_NODES[quad][j].n,\r
+ 0, v, d1, d2, d3)\r
+ - F(b * GAUSS_NODES[quad][i].n, d * GAUSS_NODES[quad][j].n,\r
+ t, 0, d1, d2, d3)\r
+ + F(b * GAUSS_NODES[quad][i].n, d * GAUSS_NODES[quad][j].n,\r
+ 0, 0, d1, d2, d3)) * GAUSS_NODES[quad][i].w\r
+ * GAUSS_NODES[quad][j].w * b * d;\r
}\r
}\r
\r
- //if (sol != sol /*|| sol == numeric_limits<double>::infinity()*/)\r
- // cout << "calcParIntQY1X1: " << sol << " isn't a Number. (" << b << ","\r
- // << d << ")(" << t << "," << v << ")(" << d1 << "," << d2 << ","\r
- // << d3 << ")" << endl;\r
-\r
return sol;\r
}\r
\r
-double calcParIntQY1(double b, double d, double t, double v, double d1,\r
- double d2, double d3) {\r
- //Fall 4\r
-\r
- double sol = 0;\r
- mexPrintf("warning: calcParIntQY1 nicht implementiert!");\r
- return 0; ///ACHTUNG noch Falsch\r
-\r
-}\r
-\r
-double calcParIntQ(double b, double d, double t, double v, double d1, double d2,\r
+/*\r
+ * Quadratur über vier Integrale\r
+ * b,d,t,v - Grenzen für Quadratur\r
+ */\r
+template<double (F)(double, double, double, double, double, double, double)>\r
+double intQ4(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
for (j = 0; j < GAUSS_SIZE[quad]; ++j) {\r
for (k = 0; k < GAUSS_SIZE[quad]; ++k) {\r
for (l = 0; l < GAUSS_SIZE[quad]; ++l) {\r
-sol += f_A(b * GAUSS_NODES[quad][i].n,\r
- d * GAUSS_NODES[quad][j].n,\r
- (t) * GAUSS_NODES[quad][k].n + d1,\r
- (v) * GAUSS_NODES[quad][l].n + d2, d3) * b\r
- * GAUSS_NODES[quad][i].w * d\r
- * GAUSS_NODES[quad][j].w * t\r
- * GAUSS_NODES[quad][k].w * v\r
- * GAUSS_NODES[quad][l].w;\r
+ sol += F(b * GAUSS_NODES[quad][i].n,\r
+ d * GAUSS_NODES[quad][j].n,\r
+ t * GAUSS_NODES[quad][k].n,\r
+ v * GAUSS_NODES[quad][l].n, d1, d2, d3) * b * d * t\r
+ * v * GAUSS_NODES[quad][i].w\r
+ * GAUSS_NODES[quad][j].w * GAUSS_NODES[quad][k].w\r
+ * GAUSS_NODES[quad][l].w;\r
+ }\r
}\r
}\r
}\r
-}\r
\r
return sol;\r
}\r
\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
-\r
-double calcOrtIntQX1X2(double b, double d, double t, double v, double d1,\r
+/*\r
+ * Grenzen in Stammfunktion über vier Integrale\r
+ * b,d,t,v - Grenzen für Stammfunktion\r
+ */\r
+template<double (F)(double, double, double, double, double, double, double)>\r
+double inline intA4(double b, double d, double t, double v, double d1,\r
double d2, double d3) {\r
- //Fall 2\r
- double sol = 0;\r
-\r
- for (int i = 0; i < GAUSS_SIZE[quad]; ++i) {\r
- for (int j = 0; j < GAUSS_SIZE[quad]; ++j) {\r
-sol += G_AY1Y2(-0.5, t + d2, v + d3, d * GAUSS_NODES[quad][j].n, 0,\r
- b * GAUSS_NODES[quad][i].n - d1) * b * d\r
- * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w;\r
- sol -= G_AY1Y2(-0.5, d2, v + d3, d * GAUSS_NODES[quad][j].n, 0,\r
- b * GAUSS_NODES[quad][i].n - d1) * b * d\r
- * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w;\r
- sol -= G_AY1Y2(-0.5, t + d2, d3, d * GAUSS_NODES[quad][j].n, 0,\r
- b * GAUSS_NODES[quad][i].n - d1) * b * d\r
- * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w;\r
- sol += G_AY1Y2(-0.5, d2, d3, d * GAUSS_NODES[quad][j].n, 0,\r
- b * GAUSS_NODES[quad][i].n - d1) * b * d\r
- * GAUSS_NODES[quad][i].w * GAUSS_NODES[quad][j].w;\r
- }\r
-}\r
\r
- return sol;\r
+ return F(b, d, t, v, d1, d2, d3) - F(b, d, t, 0, d1, d2, d3)\r
+ - F(b, d, 0, v, d1, d2, d3) + F(b, d, 0, 0, d1, d2, d3)\r
+ - F(b, 0, t, v, d1, d2, d3) + F(b, 0, t, 0, d1, d2, d3)\r
+ + F(b, 0, 0, v, d1, d2, d3) - F(b, 0, 0, 0, d1, d2, d3)\r
+ - F(0, d, t, v, d1, d2, d3) + F(0, d, t, 0, d1, d2, d3)\r
+ + F(0, d, 0, v, d1, d2, d3) - F(0, d, 0, 0, d1, d2, d3)\r
+ + F(0, 0, t, v, d1, d2, d3) - F(0, 0, t, 0, d1, d2, d3)\r
+ - F(0, 0, 0, v, d1, d2, d3) + F(0, 0, 0, 0, d1, d2, d3);\r
\r
}\r
\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
+ return intA4<F_par>(b, d, t, v, d1, d2, d3);\r
}\r
\r
double cOrtO1(double b, double d, double t, double v, double d1, double d2,\r
double d3, double* zeta) {\r
- return calcOrtIntA(b, d, t, v, d1, d2, d3);\r
+ return intA4<F_ort>(b, d, t, v, d1, d2, d3);\r
}\r
\r
double cParO2(double b, double d, double t, double v, double d1, double d2,\r
//kurze Achse nach vorn\r
// switch_dim(b, d, t, v, d1, d2);\r
\r
- if (zeta[0]*zeta[0] *(t * t + v * v) < dist2(b, d, t, v, d1, d2, d3))\r
- return calcParIntQ(b, d, t, v, d1, d2, d3);\r
+ if (zeta[0] * zeta[0] * (t * t + v * v) < dist2(b, d, t, v, d1, d2, d3))\r
+ return intQ4<f_par>(b, d, t, v, d1, d2, d3);\r
\r
- return calcParIntA(b, d, t, v, d1, d2, d3);\r
+ return intA4<F_par>(b, d, t, v, d1, d2, d3);\r
}\r
\r
double cParO3(double b, double d, double t, double v, double d1, double d2,\r
//kurze Achse nach vorn\r
// switch_dim(b, d, t, v, d1, d2);\r
\r
- if (zeta[1]*zeta[1] *(b * b + d * d) < dist2(b, d, t, v, d1, d2, d3))\r
-// if ((max(b, t) <= zeta[2] * min(d, v))\r
-// && (min(dist_s2(b, t, d1), dist_s2(d, v, d2)) > 0))\r
-// mexPrintf("Qx1x2");\r
- return calcParIntQX1X2(b, d, t, v, d1, d2, d3);\r
+ if (zeta[1] * zeta[1] * (b * b + d * d) < dist2(b, d, t, v, d1, d2, d3))\r
+ return intQ2A2<FY1Y2_par>(b, d, t, v, d1, d2, d3);\r
\r
- if (zeta[0]*zeta[0] *(t * t + v * v) < dist2(b, d, t, v, d1, d2, d3))\r
- return calcParIntQ(b, d, t, v, d1, d2, d3);\r
+ if (zeta[0] * zeta[0] * (t * t + v * v) < dist2(b, d, t, v, d1, d2, d3))\r
+ return intQ4<f_par>(b, d, t, v, d1, d2, d3);\r
\r
- return calcParIntA(b, d, t, v, d1, d2, d3);\r
+ return intA4<F_par>(b, d, t, v, d1, d2, d3);\r
}\r
\r
double cParO4(double b, double d, double t, double v, double d1, double d2,\r
switch_dim(b, d, t, v, d1, d2);\r
\r
if ((t * t + v * v) < zeta[0] * dist2(b, d, t, v, d1, d2, d3))\r
- return calcParIntQ(b, d, t, v, d1, d2, d3);\r
+ return intQ4<f_par>(b, d, t, v, d1, d2, d3);\r
\r
- if ((b * b + d * d) < zeta[0] * dist2(b, d, t, v, d1, d2, d3))\r
- return calcParIntQX1X2(b, d, t, v, d1, d2, d3);\r
+// if ((b * b + d * d) < zeta[0] * dist2(b, d, t, v, d1, d2, d3))\r
+// return calcParIntQX1X2(b, d, t, v, d1, d2, d3);\r
\r
- return calcParIntA(b, d, t, v, d1, d2, d3);\r
+ return intA4<F_par>(b, d, t, v, d1, d2, d3);\r
}\r
projects / bacc.git / commitdiff