From c78018177c2386f9dd2ad148342ecd6c163bc2ea Mon Sep 17 00:00:00 2001 From: Peter Schaefer Date: Thu, 29 Mar 2012 21:30:12 +0200 Subject: [PATCH] [src] test_v [src] Beschreibung in mex_build_V [src] zeta in Type 0 gedreht --- src/A_step.m | 21 ++++--- src/mex_build_V.cpp | 90 ++++++++++++++++++------------ src/slpRectangle.cpp | 130 +++++++++++++++++++++---------------------- src/test_v.m | 16 ++++++ 4 files changed, 146 insertions(+), 111 deletions(-) create mode 100644 src/test_v.m diff --git a/src/A_step.m b/src/A_step.m index 096c33b..0e208ce 100644 --- 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 6f11e60..27fcada 100644 --- a/src/mex_build_V.cpp +++ b/src/mex_build_V.cpp @@ -1,8 +1,26 @@ /* - * Art der Berechnung durch Parameter bestimmt... - * Analytisch/Semianalytisch/... + * MEX FUNKTION mex_build_V(coo,ele,zeta,type) * + * coo : Koordinaten Matrix + * ele : Element Matrix + * zeta : Steuerungsvariable der Zulaessigeitsbedingungen + * type : Bestimmt die Art der Berechnung * + * TYPES---- + * 1 : Voll Analytisch + * zeta wird ignoriert + * + * 0 : 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 + * + * 2 : 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 + * + * + * + * Peter Schaefer */ //#include #include @@ -10,17 +28,17 @@ #include "mex.h" /*/#include -//#include "gauss.hpp" -*/ + //#include "gauss.hpp" + */ //#define M_EPS 10e-8 /* -//#include "tbb/parallel_for.h" -*/ + //#include "tbb/parallel_for.h" + */ #include "slpRectangle.hpp" /* -//using namespace std; -//using namespace slpR; -*/ + //using namespace std; + //using namespace slpR; + */ int dimOfVec(double* vec) { if (vec[2] != 0) return 2; @@ -31,8 +49,8 @@ int dimOfVec(double* vec) { else { mexErrMsgTxt("length of Site is zero"); /* -// cerr << "dimOfVec : (" << vec[0] << " " << vec[1] << " " << vec[2] -// << ") alle Werte 0 \n";*/ + // cerr << "dimOfVec : (" << vec[0] << " " << vec[1] << " " << vec[2] + // << ") alle Werte 0 \n";*/ return -1; } @@ -44,10 +62,10 @@ inline int dimOfThird(int a, int b) { } void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { -/* -// initeQuad(); -// cout << Quad[1].nod[0]; -*/ + /* + // initeQuad(); + // cout << Quad[1].nod[0]; + */ int i, j, k; //Schleifenindizes double tmp; //Zwischenspeicherung der Einzelnen Werte int count; @@ -101,9 +119,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); // Welche Funktion soll verwendet werden - double(*ctypeP)(double, double, double, double, double, double, double, + double (*ctypeP)(double, double, double, double, double, double, double, double*); - double(*ctypeO)(double, double, double, double, double, double, double, + double (*ctypeO)(double, double, double, double, double, double, double, double*); //Art der Berechnung bestimmen @@ -261,30 +279,30 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { A[(k * em) + j] = 1. / (4 * M_PI) * tmp; if (k != j) A[(j * em) + k] = 1. / (4 * M_PI) * tmp; -/* if (count++ > ((em * (em + 1)) / 2) / 10) { - count = 0; -// mexPrintf("#"); -// cout << "#"; -// cout.flush(); - }*/ + /* if (count++ > ((em * (em + 1)) / 2) / 10) { + count = 0; + // mexPrintf("#"); + // cout << "#"; + // cout.flush(); + }*/ } } // cout << endl; //Rueckgabe (eventuell zurueck schreiben) -/* mxFree(x0); - mxFree(x1); - mxFree(x3); - mxFree(xa); - mxFree(xb); - mxFree(y0); - mxFree(y1); - mxFree(y3); - mxFree(ya); - mxFree(yb); - mxFree(d); - mxFree(dt); -*/ + /* mxFree(x0); + mxFree(x1); + mxFree(x3); + mxFree(xa); + mxFree(xb); + mxFree(y0); + mxFree(y1); + mxFree(y3); + mxFree(ya); + mxFree(yb); + mxFree(d); + mxFree(dt); + */ return; } diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp index 4937d04..6bb6781 100644 --- a/src/slpRectangle.cpp +++ b/src/slpRectangle.cpp @@ -10,7 +10,6 @@ #include "gauss.hpp" #define quad 4 // Anzahl der Quadratur Auswertungstellen 2^quad - //using namespace std; int sign(double); @@ -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)); if (sol != sol /*|| fabs(sol) == numeric_limits::infinity()*/) { - // cout << "G_AY2X2: " << sol << " isn't a Number. (" << y1 << "," << y2 - // << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl; - // cout << (x2 - y2) << "__" << sqrt(hlp) << endl; - // cout << (x1 - y1) * (x1 - y1) << " " << hlp << endl; + // cout << "G_AY2X2: " << sol << " isn't a Number. (" << y1 << "," << y2 + // << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl; + // cout << (x2 - y2) << "__" << sqrt(hlp) << endl; + // cout << (x1 - y1) * (x1 - y1) << " " << hlp << endl; } return sol; } @@ -154,8 +153,8 @@ double inline G_AY1Y2(double p, double y1, double y2, double x1, double x2, sol /= 2 * p + 2; } else { sol = NAN; - // cout << "warning in G_AY1Y2: no case for p=" << p - // << " defined. Possible: [-1.5,-0.5,0.5,...]" << endl; + // cout << "warning in G_AY1Y2: no case for p=" << p + // << " defined. Possible: [-1.5,-0.5,0.5,...]" << endl; } return sol; @@ -205,10 +204,7 @@ double inline F_ort(double x1, double x2, double y2, double y3, double d1, if ((x1 - d1) != 0) sol += (x1 - d1) - * g_AY( - 0.5, - y3, - -d3, + * g_AY(0.5, y3, -d3, sqrt( (x1 - d1) * (x1 - d1) + (x2 - y2 - d2) * (x2 - y2 - d2))); @@ -219,10 +215,7 @@ double inline F_ort(double x1, double x2, double y2, double y3, double d1, if ((y3 + d3) != 0) sol += (y3 + d3) - * g_AY( - 0.5, - x1, - d1, + * g_AY(0.5, x1, d1, sqrt( (x2 - y2 - d2) * (x2 - y2 - d2) + (y3 + d3) * (y3 + d3))); @@ -246,7 +239,7 @@ double inline F_ort(double x1, double x2, double y2, double y3, double d1, }*/ double inline apply0Int4( - double(*f)(double, double, double, double, double, double, double), + double (*f)(double, double, double, double, double, double, double), double b, double d, double t, double v, double d1, double d2, double d3) { @@ -262,7 +255,7 @@ double inline apply0Int4( } double inline apply0Int2( - double(*f)(double, double, double, double, double, double, double), + double (*f)(double, double, double, double, double, double, double), double b, double d, double t, double v, double d1, double d2, double d3) { @@ -272,8 +265,8 @@ double inline apply0Int2( } // Berechnet das eigentliche Integral fuer parallele Elemente voll Analytisch -double calcParIntA(double b, double d, double t, double v, double d1, - double d2, double d3) { +double calcParIntA(double b, double d, double t, double v, double d1, double d2, + double d3) { return apply0Int4(F_par, b, d, t, v, d1, d2, d3); } @@ -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 * gauss_nodes[quad][i].w * gauss_nodes[quad][j].w; - if (sol != sol /*|| fabs(sol) == numeric_limits::infinity()*/) { -/* 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;*/ + if (sol + != sol /*|| fabs(sol) == numeric_limits::infinity()*/) { + /* 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;*/ return sol; } @@ -380,12 +374,12 @@ double calcParIntQY1(double b, double d, double t, double v, double d1, } -double calcParIntQ(double b, double d, double t, double v, double d1, - double d2, double d3) { +double calcParIntQ(double b, double d, double t, double v, double d1, double d2, + double d3) { //Fall 0 double sol = 0; - int i,j,k,l; + int i, j, k, l; for (i = 0; i < gauss_size[quad]; ++i) { for (j = 0; j < gauss_size[quad]; ++j) { @@ -405,8 +399,8 @@ double calcParIntQ(double b, double d, double t, double v, double d1, return sol; } -double calcOrtIntA(double b, double d, double t, double v, double d1, - double d2, double d3) { +double calcOrtIntA(double b, double d, double t, double v, double d1, double d2, + double d3) { return apply0Int4(F_ort, b, d, t, v, d1, d2, d3); } @@ -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) { return calcOrtIntQX1X2(b, d, t, v, d1, d2, d3); } else { - // cout << "2error"; + // cout << "2error"; return calcOrtIntA(b, d, t, v, d1, d2, d3); } } @@ -583,8 +577,8 @@ double calcOrtIntO3(double b, double d, double t, double v, double d1, } } -double cParO0(double b, double d, double t, double v, double d1, - double d2, double d3, double* zeta) { +double cParO0(double b, double d, double t, double v, double d1, double d2, + double d3, double* zeta) { double tmp = 0; // kleine Seite nach vorn @@ -603,7 +597,7 @@ double cParO0(double b, double d, double t, double v, double d1, d2 = -d2; } - if ((t * t + v * v) * zeta[0] < d1 * d1 + d2 * d2 + d3 * d3) { + if ((t * t + v * v) < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) { return calcParIntQ(b, d, t, v, d1, d2, d3); } else { return calcParIntA(b, d, t, v, d1, d2, d3); @@ -612,13 +606,13 @@ double cParO0(double b, double d, double t, double v, double d1, return 0; } -double cParO1(double b, double d, double t, double v, double d1, - double d2, double d3, double* zeta) { +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); } -double cParO2(double b, double d, double t, double v, double d1, - double d2, double d3, double* zeta) { +double cParO2(double b, double d, double t, double v, double d1, double d2, + double d3, double* zeta) { double tmp = 0; @@ -638,7 +632,7 @@ double cParO2(double b, double d, double t, double v, double d1, d2 = -d2; } - if ((b * b + d * d) < zeta[0]* (d1 * d1 + d2 * d2 + d3 * d3)) { + if ((b * b + d * d) < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) { // cout << "E"; return calcParIntQX1X2(b, d, t, v, d1, d2, d3); } else { @@ -648,8 +642,8 @@ double cParO2(double b, double d, double t, double v, double d1, return 0; } -double cParO3(double b, double d, double t, double v, double d1, - double d2, double d3, double* zeta) { +double cParO3(double b, double d, double t, double v, double d1, double d2, + double d3, double* zeta) { double tmp = 0; @@ -669,7 +663,7 @@ double cParO3(double b, double d, double t, double v, double d1, d2 = -d2; } - if ((b * b + d * d) < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) { + if ((b * b + d * d) < zeta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) { //kleinere Achse nach vorn if (b * b + t * t > d * d + v * v) { @@ -687,7 +681,7 @@ double cParO3(double b, double d, double t, double v, double d1, d2 = tmp; } - if (min(b, t) <= zeta[1] * d1) { + if (min(b, t) <= zeta[1] * d1) { // cout << "A"; return calcParIntQY1X1(b, d, t, v, d1, d2, d3); } else { diff --git a/src/test_v.m b/src/test_v.m new file mode 100644 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 -- 2.47.3