From 5e1607ae4cb4f8eff1e3741ad1cf190b4db30601 Mon Sep 17 00:00:00 2001 From: Peter Schaefer Date: Thu, 2 Aug 2012 16:55:00 +0200 Subject: [PATCH] =?utf8?q?[src]=20F=C3=A4lle=20ummodelliert=20[src]=20Seit?= =?utf8?q?en=20Distanz=20hinzugef=C3=BCgt?= MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit --- src/mex_build_V.cpp | 15 +++++------ src/slpRectangle.cpp | 63 ++++++++++++++++++++++++-------------------- src/slpRectangle.hpp | 2 +- 3 files changed, 43 insertions(+), 37 deletions(-) diff --git a/src/mex_build_V.cpp b/src/mex_build_V.cpp index a479a0d..9d81700 100644 --- a/src/mex_build_V.cpp +++ b/src/mex_build_V.cpp @@ -11,14 +11,13 @@ * zeta wird ignoriert * * 2 : Semianalytisch im Fernfeld Quadratur ueber beide Elemente - * Quadratur bei : ((t * t + v * v) < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) : wobei v,t jeweils die laengste Seite - * in x und in y Richtung ist und d(1,2,3) der Abstandsvektor zwischen den Elementen - * - * 3 : Semianalytisch im Fernfeld Quadratur ueber ein Element - * Quadratur bei : ((b * b + d * d) < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) : wobei b,d jeweils die kuerzeste Seite - * in x und in y Richtung ist und d(1,2,3) der Abstandsvektor zwischen den Elementen + * Quadratur bei : max{dia(T1),dia(T2)} < zeta[0] * dist(T1,T2) * + * 3 : Semianalytisch wie (2) und Quadratur ueber eine Achse + * Quadratur bei : min{x1,y1} < zeta[0] * dist(x1,y1) * + * 4 : Semianalytisch wie (2) und Quadratur ueber ein Element + * Quadratur bei : min{dia(T1),dia(T2)} < zeta[0] * dist(T1,T2) * * Peter Schaefer */ @@ -135,11 +134,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { ctypeP = cParO2; ctypeO = cOrtO1; break; - case 3: // Analytisch oder Qx1x2 + case 3: // Analytisch oder Qy1x2 ctypeP = cParO3; ctypeO = cOrtO1; break; - case 4: + case 4: // Analytisch oder Qx1x2 ctypeP = cParO4; ctypeO = cOrtO1; break; diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp index d4a509a..aad61cb 100644 --- a/src/slpRectangle.cpp +++ b/src/slpRectangle.cpp @@ -73,9 +73,9 @@ void inline switch_dim(double& b, double& d, double& t, double& v, double& d1, } } -double inline dist2(double b, double d, double t, double v, double d1, double d2, - double d3) { - double dis1= 0,dis2 = 0; +double inline dist2(double b, double d, double t, double v, double d1, + double d2, double d3) { + double dis1 = 0, dis2 = 0; if (d1 < 0) dis1 = -d1 - t; else @@ -90,14 +90,28 @@ double inline dist2(double b, double d, double t, double v, double d1, double d2 if (dis2 < 0) dis2 = 0; - return dis1*dis1+dis2*dis2+d3*d3; + return dis1 * dis1 + dis2 * dis2 + d3 * d3; } - double inline dist(double b, double d, double t, double v, double d1, double d2, double d3) { return sqrt(dist2(b, d, t, v, d1, d2, d3)); } +double inline dist_s2(double b, double t, double d1) { + double dis1 = 0; + if (d1 < 0) + dis1 = -d1 - t; + else + dis1 = d1 - b; + if (dis1 < 0) + dis1 = 0; + + return dis1 * dis1; +} +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); } @@ -512,13 +526,10 @@ double cParO2(double b, double d, double t, double v, double d1, double d2, //kurze Seite nach vorn switch_site(b, d, t, v, d1, d2); - if ((t * t + v * v) < zeta[0] * dist2(b, d, t, v, d1, d2, d3)) { + if ((t * t + v * v) < zeta[0] * dist2(b, d, t, v, d1, d2, d3)) return calcParIntQ(b, d, t, v, d1, d2, d3); - } else { - return calcParIntA(b, d, t, v, d1, d2, d3); - } - return 0; + return calcParIntA(b, d, t, v, d1, d2, d3); } double cParO3(double b, double d, double t, double v, double d1, double d2, @@ -527,35 +538,31 @@ double cParO3(double b, double d, double t, double v, double d1, double d2, //kurze Seite nach vorn switch_site(b, d, t, v, d1, d2); - if ((b * b + d * d) < zeta[0] * dist2(b, d, t, v, d1, d2, d3)) { -// cout << "E"; - return calcParIntQX1X2(b, d, t, v, d1, d2, d3); - } else { - return calcParIntA(b, d, t, v, d1, d2, d3); - } + //kurze Achse nach vorn + 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); + + if ((b * b + d * d) < zeta[0] * dist_s2(b, t, d1)) + return calcParIntQY1X1(b, d, t, v, d1, d2, d3); - return 0; + return calcParIntA(b, d, t, v, d1, d2, d3); } double cParO4(double b, double d, double t, double v, double d1, double d2, double d3, double* zeta) { double tmp = 0; - double dist = dist2(b, d, t, v, d1, d2, d3); //kurze Seite nach vorn switch_site(b, d, t, v, d1, d2); - if ((t * t + v * v) < zeta[1] * dist) { - - if ((b * b + d * d) < zeta[0] * dist) { - return calcParIntQX1X2(b, d, t, v, d1, d2, d3); - } else { - return calcParIntA(b, d, t, v, d1, d2, d3); - } - } else { + 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 0; + 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); } diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp index 4995bd2..1bfd66b 100644 --- a/src/slpRectangle.hpp +++ b/src/slpRectangle.hpp @@ -45,7 +45,7 @@ double cParO1(double, double, double, double, double, double, double, double*); double cOrtO1(double, double, double, double, double, double, double, double*); //A oder Q double cParO2(double, double, double, double, double, double, double, double*); -//A oder Qx1x2 +//A oder Q oder Qy1x2 double cParO3(double, double, double, double, double, double, double, double*); // A oder Q oder Qx1x2 double cParO4(double, double, double, double, double, double, double, double*); -- 2.47.3