From: Peter Schaefer Date: Thu, 21 Mar 2013 19:37:10 +0000 (+0100) Subject: [src] dist Funktionen korrigiert X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=7d724eb02cf18a185e12ba28c90474bdec6fbeb4;p=bacc.git [src] dist Funktionen korrigiert [src] laengenparser vereinfacht --- diff --git a/src/mex_build_V.cpp b/src/mex_build_V.cpp index c361569..3e2dbfe 100644 --- a/src/mex_build_V.cpp +++ b/src/mex_build_V.cpp @@ -26,7 +26,7 @@ /***************************************************************************/ #define DEBUG 1 -#define PARALLEL +//#define PARALLEL //#include #include @@ -51,20 +51,47 @@ int inline dimOfVec(double* vec) { else if (vec[0] != 0) return 0; else { -// mexErrMsgTxt("length of Site is zero"); + mexErrMsgTxt("length of Site is zero"); return -1; } } +//// Gibt von [0 1 2] die Fehlende Zahl zurueck +//int inline dimOfThird(int a, int b) { +// return ((-(a + b) % 3) + 3) % 3; +//} + // Gibt von [0 1 2] die Fehlende Zahl zurueck int inline dimOfThird(int a, int b) { - return ((-(a + b) % 3) + 3) % 3; + switch (a + b) { + case 1: + return 2; + case 2: + return 1; + case 3: + return 0; + } } void inline add(double* a, double* b) { for (int i = 0; i < 3; ++i) a[i] += b[i]; } + +bool inline iseq(double* a, double* b) { + for (int i = 0; i < 3; ++i) + if (a[i] != b[i]) + return false; + return true; +} + +bool inline iseq(double* a, double b) { + for (int i = 0; i < 3; ++i) + if (a[i] != b) + return false; + return true; +} + void inline add(double* a, double* b, double* c) { for (int i = 0; i < 3; ++i) a[i] = b[i] + c[i]; @@ -102,19 +129,19 @@ void inline getSCorner(double* d, double x[4][3]) { } } -double inline distT(double b, double d, double t, double v, double d1, +double inline distT_par(double b, double d, double t, double v, double d1, double d2, double d3, double * dummy) { double dis1 = 0, dis2 = 0; if (d1 < 0) dis1 = -d1 - t; - else + else if(d1 > 0) dis1 = d1 - b; if (dis1 < 0) dis1 = 0; if (d2 < 0) dis2 = -d2 - v; - else + else if(d2 > 0) dis2 = d2 - d; if (dis2 < 0) dis2 = 0; @@ -123,9 +150,35 @@ double inline distT(double b, double d, double t, double v, double d1, // mexErrMsgTxt("expected d3 == 0"); // return dis1 * dis1 + dis2 * dis2 + d3 * d3; - return (dis1 * dis1 + dis2 * dis2 + d3 * d3) * (4 * M_PI); + return sqrt(dis1 * dis1 + dis2 * dis2 + d3 * d3) * (4 * M_PI); } +double inline distT_ort(double b, double d, double t, double v, double d1, + double d2, double d3, double * dummy) { + double dis1 = 0, dis2 = 0, dis3 = 0; + if (d1 < 0) + dis1 = -d1; + else + dis1 = d1 - b; + if (dis1 < 0) + dis1 = 0; + if (d2 < 0) + dis2 = -d2 - t; + else + dis2 = d2 - d; + if (dis2 < 0) + dis2 = 0; + + if (d3 < 0) + dis3 = -d3 - v; + else + dis3 = d3; + if (dis3 < 0) + dis3 = 0; + +// return dis1 * dis1 + dis2 * dis2 + dis3 * dis3; + return sqrt(dis1 * dis1 + dis2 * dis2 + dis3 * dis3) * (4 * M_PI); +} void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* mexPrintf("%d",omp_get_max_threads()); @@ -203,8 +256,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { ctypeO = cOrtO1; break; case 0: - ctypeP = distT; - ctypeO = distT; + ctypeP = distT_par; + ctypeO = distT_ort; break; } @@ -257,143 +310,165 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { #pragma omp parallel for schedule(static,1) private(i,j,k,tmp,x,y,xa,xb,ya,yb,d,dt,rx,rxa,rxb,ry,rya,ryb,lastRow,firstRow) shared(C,E,ctypeO,ctypeP,em,targetSize,actualNumberOfThreads) #endif - for (i = 0; i < actualNumberOfThreads; ++i) { - firstRow = (int) sqrt(2 * i * targetSize); - if (i == actualNumberOfThreads - 1) { - lastRow = em - 1; - } else { - lastRow = (int) sqrt(2 * (i + 1) * targetSize) - 1; - } + for (i = 0; i < actualNumberOfThreads; ++i) { + firstRow = (int) sqrt(2 * i * targetSize); + if (i == actualNumberOfThreads - 1) { + lastRow = em - 1; + } else { + lastRow = (int) sqrt(2 * (i + 1) * targetSize) - 1; + } #if DEBUG >= 4 #pragma omp critical - mexPrintf("thread Nr.: (%d) | firstRow/lastRow: [%d / %d] \n", + mexPrintf("thread Nr.: (%d) | firstRow/lastRow: [%d / %d] \n", #ifdef PARALLEL - omp_get_thread_num() + omp_get_thread_num() #else - i + i #endif - , firstRow, lastRow); + , firstRow, lastRow); #pragma omp end critical #endif - { - for (j = firstRow; j <= lastRow; ++j) { - x[0][0] = C[(int) E[j] - 1]; - x[0][1] = C[cm + (int) E[j] - 1]; - x[0][2] = C[2 * cm + (int) E[j] - 1]; + { + for (j = firstRow; j <= lastRow; ++j) { + x[0][0] = C[(int) E[j] - 1]; + x[0][1] = C[cm + (int) E[j] - 1]; + x[0][2] = C[2 * cm + (int) E[j] - 1]; - x[1][0] = C[(int) E[em + j] - 1]; - x[1][1] = C[cm + (int) E[em + j] - 1]; - x[1][2] = C[2 * cm + (int) E[em + j] - 1]; + x[1][0] = C[(int) E[em + j] - 1]; + x[1][1] = C[cm + (int) E[em + j] - 1]; + x[1][2] = C[2 * cm + (int) E[em + j] - 1]; - x[2][0] = C[(int) E[2 * em + j] - 1]; - x[2][1] = C[cm + (int) E[2 * em + j] - 1]; - x[2][2] = C[2 * cm + (int) E[2 * em + j] - 1]; + x[2][0] = C[(int) E[2 * em + j] - 1]; + x[2][1] = C[cm + (int) E[2 * em + j] - 1]; + x[2][2] = C[2 * cm + (int) E[2 * em + j] - 1]; - x[3][0] = C[(int) E[3 * em + j] - 1]; - x[3][1] = C[cm + (int) E[3 * em + j] - 1]; - x[3][2] = C[2 * cm + (int) E[3 * em + j] - 1]; + x[3][0] = C[(int) E[3 * em + j] - 1]; + x[3][1] = C[cm + (int) E[3 * em + j] - 1]; + x[3][2] = C[2 * cm + (int) E[3 * em + j] - 1]; - // Seitenvektoren aufbauen - sub(xa, x[1], x[0]); - sub(xb, x[3], x[0]); + // Seitenvektoren aufbauen + sub(xa, x[1], x[0]); + sub(xb, x[3], x[0]); - // Lageeigenschaften des Flaechenstuecks - rxa = dimOfVec(xa); - rxb = dimOfVec(xb); - rx = dimOfThird(rxa, rxb); + // Lageeigenschaften des Flaechenstuecks + rxa = dimOfVec(xa); + rxb = dimOfVec(xb); + rx = dimOfThird(rxa, rxb); + + //kleinste Ecke finden und fuer \delta verwenden +// getSCorner(dt, x); + + if (xa[rxa] <= 0 || xb[rxb] <= 0) + mexWarnMsgTxt("X Element hat Seitenlängen <=0!"); + + if(rxa == rxb) + mexWarnMsgTxt("X Laenge ist nur in einer Richtung"); + +// if (!iseq(dt, x[0])) +// mexWarnMsgTxt("x[0] ist nicht die kleinste Ecke"); + + for (k = 0; k <= j; ++k) { + + y[0][0] = C[(int) E[k] - 1]; + y[0][1] = C[cm + (int) E[k] - 1]; + y[0][2] = C[2 * cm + (int) E[k] - 1]; + + y[1][0] = C[(int) E[em + k] - 1]; + y[1][1] = C[cm + (int) E[em + k] - 1]; + y[1][2] = C[2 * cm + (int) E[em + k] - 1]; - //kleinste Ecke finden und fuer \delta verwenden - getSCorner(dt, x); + y[2][0] = C[(int) E[2 * em + k] - 1]; + y[2][1] = C[cm + (int) E[2 * em + k] - 1]; + y[2][2] = C[2 * cm + (int) E[2 * em + k] - 1]; - for (k = 0; k <= j; ++k) { + y[3][0] = C[(int) E[3 * em + k] - 1]; + y[3][1] = C[cm + (int) E[3 * em + k] - 1]; + y[3][2] = C[2 * cm + (int) E[3 * em + k] - 1]; - y[0][0] = C[(int) E[k] - 1]; - y[0][1] = C[cm + (int) E[k] - 1]; - y[0][2] = C[2 * cm + (int) E[k] - 1]; + // Seiten Vektoren aufbauen + sub(ya, y[1], y[0]); + sub(yb, y[3], y[0]); - y[1][0] = C[(int) E[em + k] - 1]; - y[1][1] = C[cm + (int) E[em + k] - 1]; - y[1][2] = C[2 * cm + (int) E[em + k] - 1]; + // Lageeigenschaften des Flaechenstuecks + rya = dimOfVec(ya); + ryb = dimOfVec(yb); + ry = dimOfThird(rya, ryb); + + //kleinste Ecke finden +// getSCorner(d, y); - y[2][0] = C[(int) E[2 * em + k] - 1]; - y[2][1] = C[cm + (int) E[2 * em + k] - 1]; - y[2][2] = C[2 * cm + (int) E[2 * em + k] - 1]; + if (ya[rya] <= 0 || yb[ryb] <= 0) + mexWarnMsgTxt("Y Element hat Seitenlängen <=0!"); - y[3][0] = C[(int) E[3 * em + k] - 1]; - y[3][1] = C[cm + (int) E[3 * em + k] - 1]; - y[3][2] = C[2 * cm + (int) E[3 * em + k] - 1]; + if(rya == ryb) + mexWarnMsgTxt("Y Laenge ist nur in einer Richtung"); - // Seiten Vektoren aufbauen - sub(ya, y[1], y[0]); - sub(yb, y[3], y[0]); +// if (!iseq(d, y[0])) +// mexWarnMsgTxt("y[0] ist nicht die kleinste Ecke"); - // Lageeigenschaften des Flaechenstuecks - rya = dimOfVec(ya); - ryb = dimOfVec(yb); - ry = dimOfThird(rya, ryb); + //\delta berechnen +// sub(d, dt); - //kleinste Ecke finden - getSCorner(d, y); + sub(d, y[0], x[0]); + + if(j==1 && k==0){ + mexPrintf("[%d %d](%f,%f)(%f,%f)+(%f,%f,%f)", j, k, + (xa[rxa]), (xb[rxb]), (ya[rya]), + (yb[ryb]), d[rxa], d[rxb], d[rx]); + mexPrintf("{%d,%d;%d,%d}\n",rxa,rxb,rya,rxb); + } - //\delta berechnen - sub(d, dt); #if DEBUG >= 5 #pragma omp critical - mexPrintf("[%d %d](%f,%f)(%f,%f)+(%f,%f,%f)\n", j, k, - fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]), - fabs(yb[rxb]), d[rxa], d[rxb], d[rx]); + mexPrintf("[%d %d](%f,%f)(%f,%f)+(%f,%f,%f)\n", j, k, + (xa[rxa]), (xb[rxb]), (ya[rya]), + (yb[ryb]), d[rxa], d[rxb], d[rx]); #pragma omp end critical #endif - if (rx == ry) { - if (rxa == rya) { - tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), - fabs(ya[rxa]), fabs(yb[rxb]), d[rxa], - d[rxb], d[rx], zeta); - - } else { - tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), - fabs(yb[rxa]), fabs(ya[rxb]), d[rxa], - d[rxb], d[rx], zeta); - } + if (rx == ry) { + if (rxa == rya) { + tmp = ctypeP((xa[rxa]), (xb[rxb]), (ya[rxa]), + (yb[rxb]), d[rxa], d[rxb], d[rx], zeta); } else { + tmp = ctypeP((xa[rxa]), (xb[rxb]), (yb[rxa]), + (ya[rxb]), d[rxa], d[rxb], d[rx], zeta); + } - if (rxa == rya) { - tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), - fabs(ya[rya]), fabs(yb[ryb]), d[rxb], - d[rxa], d[rx], zeta); - } else if (rxa == ryb) { - tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), - fabs(yb[ryb]), fabs(ya[rya]), d[rxb], - d[rxa], d[rx], zeta); - } else if (rxb == rya) { - tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), - fabs(ya[rya]), fabs(yb[ryb]), d[rxa], - d[rxb], d[rx], zeta); - } else { - tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), - fabs(yb[ryb]), fabs(ya[rya]), d[rxa], - d[rxb], d[rx], zeta); - } - + } else { + + if (rxa == rya) { + tmp = ctypeO((xb[rxb]), (xa[rxa]), (ya[rya]), + (yb[ryb]), d[rxb], d[rxa], d[rx], zeta); + } else if (rxa == ryb) { + tmp = ctypeO((xb[rxb]), (xa[rxa]), (yb[ryb]), + (ya[rya]), d[rxb], d[rxa], d[rx], zeta); + } else if (rxb == rya) { + tmp = ctypeO((xa[rxa]), (xb[rxb]), (ya[rya]), + (yb[ryb]), d[rxa], d[rxb], d[rx], zeta); + } else { + tmp = ctypeO((xa[rxa]), (xb[rxb]), (yb[ryb]), + (ya[rya]), d[rxa], d[rxb], d[rx], zeta); } + } + #pragma omp critical - A[(k * em) + j] = 1. / (4 * M_PI) * tmp; - if (k != j) - A[(j * em) + k] = 1. / (4 * M_PI) * tmp; + A[(k * em) + j] = 1. / (4 * M_PI) * tmp; + if (k != j) + A[(j * em) + k] = 1. / (4 * M_PI) * tmp; #pragma omp end critical - } } } } -#ifdef PARALLEL } +#ifdef PARALLEL +} #endif return; diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp index 9828f44..04a0c8a 100644 --- a/src/slpRectangle.cpp +++ b/src/slpRectangle.cpp @@ -97,7 +97,7 @@ void inline switch_dim(double& b, double& d, double& t, double& v, double& d1, /* * gibt den Abstand zum Quadrat zurueck */ -double inline dist2(double b, double d, double t, double v, double d1, +double inline dist2_par(double b, double d, double t, double v, double d1, double d2, double d3) { double dis1 = 0, dis2 = 0; if (d1 < 0) @@ -116,14 +116,38 @@ double inline dist2(double b, double d, double t, double v, double d1, return dis1 * dis1 + dis2 * dis2 + d3 * d3; } + /* - * gibt den Abstand zurueck + * gibt den Abstand zum Quadrat zurueck */ -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 dist2_ort(double b, double d, double t, double v, double d1, + double d2, double d3) { + double dis1 = 0, dis2 = 0, dis3 = 0; + if (d1 < 0) + dis1 = -d1; + else + dis1 = d1 - b; + if (dis1 < 0) + dis1 = 0; + + if (d2 < 0) + dis2 = -d2 - t; + else + dis2 = d2 - d; + if (dis2 < 0) + dis2 = 0; + + if (d3 < 0) + dis3 = -d3 - v; + else + dis3 = d3; + if (dis3 < 0) + dis3 = 0; + + return dis1 * dis1 + dis2 * dis2 + dis3 * dis3; } + /* * gibt den Abstand in einer Richtung zum Quadrat zurueck */ @@ -138,12 +162,7 @@ double inline dist_s2(double b, double t, double d1) { return dis1 * dis1; } -/* - * gibt den Abstand in einer Richtung zurueck - */ -double inline dist_s(double b, double t, double d1) { - return sqrt(dist_s2(b, t, d1)); -} + double inline f_par(double x1, double x2, double y1, double y2, double d1, double d2, double d3) { @@ -288,7 +307,7 @@ double inline G_AY1Y2(double p, double y1, double y2, double x1, double x2, double inline Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2, double l) { double sol = 0; - mexPrintf("warning: Gs_AX2Y1Y2 nicht implementiert!"); + mexWarnMsgTxt("Gs_AX2Y1Y2 nicht implementiert!"); return sol; } @@ -320,7 +339,7 @@ double inline F_par(double x1, double x2, double y1, double y2, double d1, double inline F_ort(double x1, double x2, double y2, double y3, double d1, double d2, double d3) { - // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y2,y3,d1,d2,d3); +// mexPrintf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f \n",x1,x2,y2,y3,d1,d2,d3); double sol = -G_AY1Y2(0.5, y3, x1, -d3, d1, x2 - y2 - d2); if ((x1 - d1) * (x2 - y2 - d2) != 0) @@ -348,19 +367,16 @@ double inline F_ort(double x1, double x2, double y2, double y3, double d1, return sol / 2.; } - -double inline FY1Y2_par(double x1, double x2, double y1, double y2, double d1, double d2, - double d3) { +double inline FY1Y2_par(double x1, double x2, double y1, double y2, double d1, + double d2, double d3) { return G_AY1Y2(-0.5, y1 + d1, y2 + d2, x1, x2, d3); } -double inline FY2X2_par(double y1, double x1, double y2, double x2, double d1, double d2, - double d3) { +double inline FY2X2_par(double y1, double x1, double y2, double x2, double d1, + double d2, double d3) { return G_AY2X2(y1 + d1, y2 + d2, x1, x2, d3); } - - /* * Quadratur über zwei Integrale und Grenzen in Stammfunktion über zwei Integrale * b,d - Grenzen für Quadratur @@ -439,11 +455,15 @@ double inline intA4(double b, double d, double t, double v, double d1, double cParO1(double b, double d, double t, double v, double d1, double d2, double d3, double* zeta) { +// if (d3 && fabs(d3)!=2) +// mexPrintf("%f \n", 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) { +// if (!d3) +// mexPrintf("%f \n", d3); return intA4(b, d, t, v, d1, d2, d3); } @@ -457,13 +477,12 @@ 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)) + if (zeta[0] * zeta[0] * (t * t + v * v) < dist2_par(b, d, t, v, d1, d2, d3)) return intQ4(b, d, t, v, d1, d2, d3); return intA4(b, d, t, v, d1, d2, d3); } - double cOrtO2(double b, double d, double t, double v, double d1, double d2, double d3, double* zeta) { double tmp = 0, dis2 = 0; @@ -481,23 +500,15 @@ double cOrtO2(double b, double d, double t, double v, double d1, double d2, tmp = b; b = v; v = tmp; - d1 = -d1; - d3 = -d3; + tmp = -d1; + d1 = -d3; + d3 = tmp; } - if (d2 < 0) - tmp = -d2 - t; - else - tmp = d2 - d; - if (tmp < 0) - tmp = 0; - - dis2 = tmp * tmp + d1 * d1 + d3 * d3; - //kurze Achse nach vorn // switch_dim(b, d, t, v, d1, d2); - if (zeta[0] * zeta[0] * (t * t + v * v) < dis2) + if (zeta[0] * zeta[0] * (t * t + v * v) < dist2_ort(b,d,t,v,d1,d2,d3)) return intQ4(b, d, t, v, d1, d2, d3); return intA4(b, d, t, v, d1, d2, d3); @@ -512,10 +523,10 @@ 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 (zeta[1] * zeta[1] * (b * b + d * d) < dist2_par(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)) + if (zeta[0] * zeta[0] * (t * t + v * v) < dist2_par(b, d, t, v, d1, d2, d3)) return intQ4(b, d, t, v, d1, d2, d3); return intA4(b, d, t, v, d1, d2, d3); @@ -532,10 +543,10 @@ double cParO4(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 ((t * t + v * v) < zeta[0] * dist2(b, d, t, v, d1, d2, d3)) + if ((t * t + v * v) < zeta[0] * dist2_par(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)) +// if ((b * b + d * d) < zeta[0] * dist2_par(b, d, t, v, d1, d2, d3)) // return calcParIntQX1X2(b, d, t, v, d1, d2, d3); return intA4(b, d, t, v, d1, d2, d3);