From: treecity Date: Wed, 18 Jan 2012 10:03:26 +0000 (+0000) Subject: [src] A_plots automatische Plot Generierung X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=1e3a6b7c77f48934ff570f34c2302ebe45950b32;p=bacc.git [src] A_plots automatische Plot Generierung [src] mex Skript eta und eps angepasst git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@74 26120e32-c555-405d-b3e1-1f783fb42516 --- diff --git a/.cproject b/.cproject index b3be917..dbf0e66 100644 --- a/.cproject +++ b/.cproject @@ -12,11 +12,12 @@ - - - + + + + @@ -31,41 +32,34 @@ - + - + - - + - + + - - + - @@ -78,7 +72,15 @@ - + + + + + + + + + @@ -96,49 +98,26 @@ - + - + - + + - - + - diff --git a/.settings/org.eclipse.ltk.core.refactoring.prefs b/.settings/org.eclipse.ltk.core.refactoring.prefs new file mode 100644 index 0000000..044f491 --- /dev/null +++ b/.settings/org.eclipse.ltk.core.refactoring.prefs @@ -0,0 +1,3 @@ +#Mon Jan 16 22:45:22 CET 2012 +eclipse.preferences.version=1 +org.eclipse.ltk.core.refactoring.enable.project.refactoring.history=false diff --git a/src/A_loadMeshF.m b/src/A_loadMeshF.m index fa59a80..9879040 100644 --- a/src/A_loadMeshF.m +++ b/src/A_loadMeshF.m @@ -21,7 +21,7 @@ function A_loadMeshF(file) if(exist('time','var')~=0) G_T = time; end - if(exist('datA','var')~=0) + if(exist('data','var')~=0) G_D = data; end end diff --git a/src/A_plot.m b/src/A_plot.m index 72a1dc0..dffc3a9 100644 --- a/src/A_plot.m +++ b/src/A_plot.m @@ -15,16 +15,18 @@ if step<1 else figure(4) -loglog(G_D(:,1),G_D(:,[3+(0:step-1)*3])) +loglog(G_D(:,1),[G_D(:,[3+(0:step-1)*3]),20*G_D(:,1).^(-1/2),6*G_D(:,1).^(-1/4),5*G_D(:,1).^(-3/4)],... + G_D(1:end-1,1), sqrt(G_D(2:end,[4+(0:step-1)*3])-G_D(1:end-1,[4+(0:step-1)*3])),'--') +title('Fehler') xlabel('Elemente'); ylabel('Schaetzer'); -legend(type2str (G_D(1,[2+(0:step-1)*3])',:)); +legend([['m';'m';'m'] type2str(G_D(1,[2+(0:step-1)*3])',:); 'N^{-1/2} '; 'N^{-1/4} '; 'N^{-3/4} ';['g';'g';'g'] type2str(G_D(1,[2+(0:step-1)*3])',:)] ,'location','best'); figure(5) loglog(G_D(:,1),G_D(:,[4+(0:step-1)*3])) - +title('Energie Norm') xlabel('Elemente'); ylabel('eNorm'); legend(type2str(G_D(1,[2+(0:step-1)*3])',:)); diff --git a/src/A_plots.m b/src/A_plots.m new file mode 100644 index 0000000..deecfe7 --- /dev/null +++ b/src/A_plots.m @@ -0,0 +1,90 @@ +function A_plots(files,printt) + +type2str = {'Analytisch ' 'Quad Element' 'Quad Achse ' 'Quad Seite '}; + +G_D=[]; +X = []; +leg1 = {}; +leg2 = {}; +leg3 = {}; + +for i = 1:length(files) + + load(files{i}); + + [m n] = size(data); + step = round(n/3); + if step<1 + disp (['Error: No Data to show. : ' i]) + continue; + end + + if(size(X,1)>m&&size(X,1)~=0) + data(size(X,1),end) = 0; + end + if(size(X,1) - -using namespace std; - -typedef struct _Gauss{ - double n; - double w; -} Gauss; - -Gauss gauss[][32] = { - { - {0.5,1}}, - { - {0.2113248654051871,0.4999999999999999}, - {0.7886751345948129,0.4999999999999999}}, - { - {0.06943184420297371,0.1739274225687271}, - {0.3300094782075718,0.3260725774312732}, - {0.6699905217924281,0.326072577431273}, - {0.9305681557970263,0.173927422568727}}, - { - {0.01985507175123191,0.05061426814518822}, - {0.1016667612931866,0.1111905172266873}, - {0.2372337950418356,0.1568533229389437}, - {0.4082826787521751,0.1813418916891811}, - {0.5917173212478251,0.1813418916891812}, - {0.7627662049581644,0.1568533229389436}, - {0.8983332387068135,0.1111905172266873}, - {0.980144928248768,0.05061426814518825}}, - { - {0.005299532504175253,0.01357622970587728}, - {0.02771248846338364,0.03112676196932351}, - {0.06718439880608407,0.04757925584124631}, - {0.1222977958224983,0.06231448562776711}, - {0.1910618777986778,0.07479799440828819}, - {0.2709916111713863,0.08457825969750139}, - {0.3591982246103706,0.09130170752246231}, - {0.4524937450811814,0.09472530522753414}, - {0.5475062549188187,0.09472530522753445}, - {0.6408017753896295,0.09130170752246193}, - {0.7290083888286135,0.08457825969750142}, - {0.808938122201322,0.0747979944082887}, - {0.8777022041775016,0.06231448562776742}, - {0.9328156011939159,0.04757925584124623}, - {0.9722875115366163,0.03112676196932403}, - {0.9947004674958249,0.0135762297058772}}, - { - {0.001368069075259104,0.003509305004735008}, - {0.007194244227365809,0.008137197365452637}, - {0.01761887220624681,0.01269603265463125}, - {0.03254696203113011,0.01713693145651049}, - {0.05183942211697379,0.02141794901111367}, - {0.07531619313371507,0.02549902963118811}, - {0.1027581020160289,0.0293420467392677}, - {0.1339089406298552,0.03291111138818078}, - {0.1684778665348924,0.03617289705442425}, - {0.2061421213796186,0.03909694789353495}, - {0.2465500455338854,0.04165596211347347}, - {0.2893243619346824,0.04382604650220204}, - {0.3340656988589361,0.04558693934788215}, - {0.3803563188739316,0.04692219954040191}, - {0.4277640192086017,0.04781936003963715}, - {0.4758461671561307,0.04827004425736407}, - {0.5241538328438691,0.04827004425736345}, - {0.5722359807913984,0.04781936003963707}, - {0.6196436811260685,0.04692219954040222}, - {0.6659343011410638,0.04558693934788176}, - {0.7106756380653179,0.04382604650220195}, - {0.7534499544661147,0.04165596211347328}, - {0.7938578786203814,0.03909694789353577}, - {0.8315221334651075,0.03617289705442422}, - {0.866091059370145,0.03291111138818064}, - {0.8972418979839712,0.02934204673926768}, - {0.9246838068662852,0.02549902963118814}, - {0.9481605778830262,0.0214179490111134}, - {0.9674530379688698,0.01713693145651111}, - {0.9823811277937532,0.01269603265463138}, - {0.992805755772634,0.008137197365453092}, - {0.9986319309247407,0.003509305004734969}}}; - -int main(){ - - cout << "--" << endl; - cout << gauss[0][0].w << endl; - return 0; -} - -#endif diff --git a/src/gauss.hpp b/src/gauss.hpp new file mode 100644 index 0000000..0d68bbf --- /dev/null +++ b/src/gauss.hpp @@ -0,0 +1,121 @@ +#ifndef GAUSS_NODES +#define GAUSS_NODES +/* + * Gaus Auswertungspunkte auf dem Intervall [0 1] + * + * 2^$NUM entspricht der Anzahl der Auswertungspunkte + * + * Anzahl der Punkte: gauss_size[$NUM] + * Auswertungstellen: gauss_nodes[$NUM][$I].n ($I aktuelle Position) + * Gewichte: gauss_nodes[$NUM][$I].w ($I aktuelle Position) + * + * Beispiel: + * + * > double sol; + * > for(int i=0;i sol += sin(gauss_nodes[3][i].n) * gauss_nodes[3][i].w; + * + * Peter Schaefer + */ + +typedef struct _gauss{ + double n; + double w; +} gauss; + +/* +int inline exp2(int exp){ + int e = 1; + for(int i=0; i + +using namespace std; + +int main(){ + + cout << "--" << endl; + cout << gauss[0][0].w << endl; + return 0; +} +*/ +#endif diff --git a/src/mark.m b/src/mark.m index aae5d94..83af46f 100644 --- a/src/mark.m +++ b/src/mark.m @@ -1,5 +1,5 @@ -function REF = mark(xF2S,ind,eta,eps) -% function REF = mark(xF2S,ind,eta,eps) +function REF = mark(xF2S,ind,theta,eta) +% function REF = mark(xF2S,ind,theta,eta) % xF2S - Father son % ind - error estimator % eta - refine element? (0..1, 1 = All) @@ -18,12 +18,12 @@ t3 =0; t4 = 0; Ct = T4*xF2S; %% Muss ueberhaupt verfeinert werden -if(eta <1) +if(theta <1) [s_ind idx] = sort(ind,'descend'); sum_ind = cumsum(s_ind,1); - ell = find(sum_ind >= sum_ind(end) * eta,1); + ell = find(sum_ind >= sum_ind(end) * theta,1); t1(idx(ell:end)) = 1; end @@ -31,9 +31,9 @@ end %% Wie muss verfeinert werden if(eps > 0) - t3 = (eps*abs(Ct(3,:)) >= sqrt(Ct(2,:).^2+Ct(4,:).^2)); + t3 = (eta*abs(Ct(3,:)) >= sqrt(Ct(2,:).^2+Ct(4,:).^2)); REF(t3) = 3; - t4 = (eps*abs(Ct(4,:)) >= sqrt(Ct(2,:).^2+Ct(3,:).^2)); + t4 = (eta*abs(Ct(4,:)) >= sqrt(Ct(2,:).^2+Ct(3,:).^2)); REF(t4) = 4; end REF(~(t4+t3+t1)) = 2; diff --git a/src/marked.m b/src/marked.m new file mode 100644 index 0000000..e69de29 diff --git a/src/mesh.tar.gz b/src/mesh.tar.gz new file mode 100644 index 0000000..0a432f9 Binary files /dev/null and b/src/mesh.tar.gz differ diff --git a/src/mex_build_AU.cpp b/src/mex_build_AU.cpp index 7f54061..2648cd0 100644 --- a/src/mex_build_AU.cpp +++ b/src/mex_build_AU.cpp @@ -17,21 +17,21 @@ //#include "tbb/parallel_for.h" - #include "slpRectangle.hpp" using namespace std; using namespace slpR; int dimOfVec(double* vec) { + std::cout << ""; if (vec[2] != 0) return 2; - else if (vec[1] !=0) + else if (vec[1] != 0) return 1; - else if (vec[0] !=0) + else if (vec[0] != 0) return 0; else { - cerr << "dimOfVec : (" << vec[0] << " " << vec[1]<< " " << vec[2] + cerr << "dimOfVec : (" << vec[0] << " " << vec[1] << " " << vec[2] << ") alle Werte 0 \n"; return -1; } @@ -77,7 +77,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double * C = mxGetPr(prhs[0]); double * E = mxGetPr(prhs[1]); - int type = (int) *(mxGetPr(prhs[3])); double * mu = mxGetPr(prhs[2]); @@ -113,6 +112,10 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { ctypeP = cParO1; ctypeO = NULL; break; + case 0: + ctypeP = cParO0; + ctypeO = NULL; + break; case 2: ctypeP = cParO2; ctypeO = NULL; @@ -178,7 +181,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { } } - for (k = 0; k < em; ++k) { + for (k = 0; k <= j; ++k) { y0[0] = C[(int) E[k] - 1]; y0[1] = C[cm + (int) E[k] - 1]; y0[2] = C[2 * cm + (int) E[k] - 1]; @@ -229,50 +232,39 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 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], mu); + tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]), + fabs(yb[rxb]), d[rxa], d[rxb], d[rx], mu); } else { - tmp - = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), - fabs(yb[rxa]), fabs(ya[rxb]), d[rxa], - d[rxb], d[rx], mu); + tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[rxa]), + fabs(ya[rxb]), d[rxa], d[rxb], d[rx], mu); } } else { if (rxa == rya) { - tmp - = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), - fabs(ya[rya]), fabs(yb[ryb]), d[rxb], - d[rxa], d[rx], mu); + tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), fabs(ya[rya]), + fabs(yb[ryb]), d[rxb], d[rxa], d[rx], mu); } else if (rxa == ryb) { - tmp - = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), - fabs(yb[ryb]), fabs(ya[rya]), d[rxb], - d[rxa], d[rx], mu); + tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), fabs(yb[ryb]), + fabs(ya[rya]), d[rxb], d[rxa], d[rx], mu); } else if (rxb == rya) { - tmp - = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), - fabs(ya[rya]), fabs(yb[ryb]), d[rxa], - d[rxb], d[rx], mu); + tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rya]), + fabs(yb[ryb]), d[rxa], d[rxb], d[rx], mu); } else { - tmp - = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), - fabs(yb[ryb]), fabs(ya[rya]), d[rxa], - d[rxb], d[rx], mu); + tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[ryb]), + fabs(ya[rya]), d[rxa], d[rxb], d[rx], mu); } } 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){ + if (k != j) + A[(j * em) + k] = 1. / (4 * M_PI) * tmp; +/* if (count++ > ((em * (em + 1)) / 2) / 10) { count = 0; - cout << "#"; - cout.flush(); +// mexPrintf("#"); +// cout << "#"; +// cout.flush(); }*/ } @@ -280,14 +272,18 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { } cout << endl; //Rueckgabe (eventuell zurueck schreiben) - //mxFree(x0); - //mxFree(x1); - //mxFree(x3); - //mxFree(y0); - //mxFree(y1); - //mxFree(y3); - //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 ccdfc86..a6d212c 100644 --- a/src/slpRectangle.cpp +++ b/src/slpRectangle.cpp @@ -7,21 +7,20 @@ //#include #include "slpRectangle.hpp" -//#include "gauss.hpp" +#include "gauss.hpp" +#define quad 2 // Anzahl der Quadratur Auswertungstellen 2^quad using namespace std; //using namespace slpR; -double gauss_w5[] = { 0.1185, 0.2393, 0.2844, 0.2393, 0.1185 }; -double gauss_n5[] = { 0.0469, 0.2308, 0.5000, 0.7692, 0.9531 }; +//double gauss_w5[] = { 0.1185, 0.2393, 0.2844, 0.2393, 0.1185 }; +//double gauss_n5[] = { 0.0469, 0.2308, 0.5000, 0.7692, 0.9531 }; int sign(double); //double arsinh(double); - /* ============================== */ - int inline sign(double x) { return x > 0 ? 1 : (x < 0 ? -1 : 0); } @@ -39,13 +38,19 @@ double inline min(double x, double y) { return x > y ? y : x; } +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 slpR::g_QY(double p, double y, double x, double l) { double sol = 0; - for (int i = 0; i < 5; ++i) { - sol += pow((x - gauss_n5[i] * y) * (x - gauss_n5[i] * y) + l * l, p) - * gauss_w5[i] * y; + 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; } return sol; @@ -59,8 +64,8 @@ double slpR::g_AY(double p, double y, double x, double l) { if (l != 0) { if (p == 0.5) { - sol = (y - x) / 2 * sqrt((y - x) * (y - x) + l * l) + l * l / 2 - * asinh((y - x) / fabs(l)); + sol = (y - x) / 2 * sqrt((y - x) * (y - x) + l * l) + + l * l / 2 * asinh((y - x) / fabs(l)); // printf("%.2f |",sol); } else if (p == 0) sol = y - x; @@ -71,8 +76,8 @@ double slpR::g_AY(double p, double y, double x, double l) { else if (p == -1.5) sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l)); else - sol = (y - x) * pow((y - x) * (y - x) + l * l, p) + 2 * p * l * l - * g_AY(p - 1, y, x, l) / (2 * p + 1); + sol = (y - x) * pow((y - x) * (y - x) + l * l, p) + + 2 * p * l * l * g_AY(p - 1, y, x, l) / (2 * p + 1); } else { if (p == -0.5) sol = sign(y - x) * log(fabs(y - x)); @@ -83,15 +88,19 @@ double slpR::g_AY(double p, double y, double x, double l) { return sol; } -double slpR::G_QY1Y2(double p, double y1, double y2, double x1, double x2, double l) { +double slpR::G_QY1Y2(double p, double y1, double y2, double x1, double x2, + double l) { double sol = 0; - for (int i = 0; i < 5; ++i) { - for (int j = 0; j < 5; ++j) { + for (int i = 0; i < gauss_size[quad]; ++i) { + for (int j = 0; j < gauss_size[quad]; ++j) { sol += pow( - (x1 - y1 * gauss_n5[i]) * (x1 - y1 * gauss_n5[i]) + (x2 - - y2 * gauss_n5[j]) * (x2 - y2 * gauss_n5[j]) + l - * l, p) * y1 * gauss_w5[i] * y2 * gauss_w5[j]; + (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; @@ -105,44 +114,46 @@ double slpR::G_AY2X2(double y1, double y2, double x1, double x2, double d3) { double sol = sqrt(hlp); - if ((x2 - y2) != 0) //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)==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; + 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; } return sol; } -double slpR::G_AY1Y2(double p, double y1, double y2, double x1, double x2, double l) { +double slpR::G_AY1Y2(double p, double y1, double y2, double x1, double x2, + double l) { // printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l); double pt = p + 1.5; double sol = 0; if (pt == 0) { if (l == 0) { - sol = -sqrt((y1 - x1) * (y1 - x1) + (y2 - x2) * (y2 - x2)) / ((y1 - - x1) * (y2 - x2)); + sol = -sqrt((y1 - x1) * (y1 - x1) + (y2 - x2) * (y2 - x2)) + / ((y1 - x1) * (y2 - x2)); } else { - sol = sign((y1 - x1) * (y2 - x2)) / (2 * fabs(l)) * acos( - -2 * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2) / (((y1 - - x1) * (y1 - x1) + l * l) * ((y2 - x2) * (y2 - x2) - + l * l)) + 1); + sol = sign((y1 - x1) * (y2 - x2)) / (2 * fabs(l)) + * acos( + -2 * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2) + / (((y1 - x1) * (y1 - x1) + l * l) + * ((y2 - x2) * (y2 - x2) + l * l)) + + 1); } } else if ((pt > 0) && ((int) pt == pt)) { if (l != 0) sol = 2 * p * l * l * G_AY1Y2(p - 1, y1, y2, x1, x2, l); if ((y1 - x1) != 0) - sol += (y1 - x1) * g_AY(p, y2, x2, - sqrt((y1 - x1) * (y1 - x1) + l * l)); + sol += (y1 - x1) + * g_AY(p, y2, x2, sqrt((y1 - x1) * (y1 - x1) + l * l)); if ((y2 - x2) != 0) - sol += (y2 - x2) * g_AY(p, y1, x1, - sqrt((y2 - x2) * (y2 - x2) + l * l)); + sol += (y2 - x2) + * g_AY(p, y1, x1, sqrt((y2 - x2) * (y2 - x2) + l * l)); sol /= 2 * p + 2; } else { sol = NAN; @@ -160,8 +171,8 @@ double slpR::Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2, return sol; } -double slpR::F_par(double x1, double x2, double y1, double y2, double d1, double d2, - double d3) { +double slpR::F_par(double x1, double x2, double y1, double y2, double d1, + double d2, double d3) { // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3); double sol = (x1 - y1 - d1) * (x2 - y2 - d2); @@ -170,40 +181,54 @@ double slpR::F_par(double x1, double x2, double y1, double y2, double d1, double sol *= G_AY1Y2(-0.5, x1, x2, y1 + d1, y2 + d2, d3); if ((x1 - y1 - d1) != 0) - sol -= (x1 - y1 - d1) * g_AY(0.5, x1, y1 + d1, - sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3)); + sol -= (x1 - y1 - d1) + * g_AY(0.5, x1, y1 + d1, + sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3)); if ((x2 - y2 - d2) != 0) - sol -= (x2 - y2 - d2) * g_AY(0.5, x2, y2 + d2, - sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3)); + sol -= (x2 - y2 - d2) + * g_AY(0.5, x2, y2 + d2, + sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3)); - double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2 - - d2) + d3 * d3); + double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + + (x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3); sol += 1. / 3 * hlp * sqrt(hlp); return sol; } -double slpR::F_ort(double x1, double x2, double y2, double y3, double d1, double d2, - double d3) { +double slpR::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); double sol = -G_AY1Y2(0.5, y3, x1, -d3, d1, x2 - y2 - d2); if ((x1 - d1) * (x2 - y2 - d2) != 0) - sol -= (x1 - d1) * (x2 - y2 - d2) * G_AY1Y2(-0.5, x2, y3, y2 + d2, -d3, - x1 - d1); + sol -= (x1 - d1) * (x2 - y2 - d2) + * G_AY1Y2(-0.5, x2, y3, y2 + d2, -d3, x1 - d1); if ((x1 - d1) != 0) - sol += (x1 - d1) * g_AY(0.5, y3, -d3, - sqrt((x1 - d1) * (x1 - d1) + (x2 - y2 - d2) * (x2 - y2 - d2))); + sol += (x1 - d1) + * g_AY( + 0.5, + y3, + -d3, + sqrt( + (x1 - d1) * (x1 - d1) + + (x2 - y2 - d2) * (x2 - y2 - d2))); if ((y3 + d3) * (x2 - y2 - d2) != 0) - sol -= (y3 + d3) * (x2 - y2 - d2) * G_AY1Y2(-0.5, x1, x2, d1, y2 + d2, - -y3 - d3); + sol -= (y3 + d3) * (x2 - y2 - d2) + * G_AY1Y2(-0.5, x1, x2, d1, y2 + d2, -y3 - d3); if ((y3 + d3) != 0) - sol += (y3 + d3) * g_AY(0.5, x1, d1, - sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + (y3 + d3) * (y3 + d3))); + sol += (y3 + d3) + * g_AY( + 0.5, + x1, + d1, + sqrt( + (x2 - y2 - d2) * (x2 - y2 - d2) + + (y3 + d3) * (y3 + d3))); return sol / 2.; } @@ -225,24 +250,27 @@ double slpR::F_ort(double x1, double x2, double y2, double y3, double d1, double double slpR::apply0Int4( double(*f)(double, double, double, double, double, double, double), - double b, double d, double t, double v, double d1, double d2, double d3) { + 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); + 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 slpR::apply0Int2( double(*f)(double, double, double, double, double, double, double), - double b, double d, double t, double v, double d1, double d2, double d3) { + 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, 0, t, v, d1, d2, d3) - f(0, d, t, - v, d1, d2, d3) + f(0, 0, t, v, d1, d2, 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); } @@ -258,16 +286,20 @@ double slpR::calcParIntQX1X2(double b, double d, double t, double v, double d1, //Fall 2 double sol = 0; - for (int i = 0; i < 5; ++i) { - for (int j = 0; j < 5; ++j) { - sol += G_AY1Y2(-0.5, t + d1, v + d2, b * gauss_n5[i], - d * gauss_n5[j], d3) * b * d * gauss_w5[i] * gauss_w5[j]; - sol -= G_AY1Y2(-0.5, d1, v + d2, b * gauss_n5[i], d * gauss_n5[j], - d3) * b * d * gauss_w5[i] * gauss_w5[j]; - sol -= G_AY1Y2(-0.5, t + d1, d2, b * gauss_n5[i], d * gauss_n5[j], - d3) * b * d * gauss_w5[i] * gauss_w5[j]; - sol += G_AY1Y2(-0.5, d1, d2, b * gauss_n5[i], d * gauss_n5[j], d3) - * b * d * gauss_w5[i] * gauss_w5[j]; + 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; } } @@ -280,40 +312,51 @@ double slpR::calcParIntQY1X1(double b, double d, double t, double v, double d1, //Fall 3 double sol = 0; - for (int i = 0; i < 5; ++i) { - for (int j = 0; j < 5; ++j) { - - - sol - += G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], d, d3) - * t * b * gauss_w5[i] * gauss_w5[j]; - sol - -= G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], 0, d3) - * t * b * gauss_w5[i] * gauss_w5[j]; - sol -= G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], d, d3) - * t * b * gauss_w5[i] * gauss_w5[j]; - sol += G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], 0, d3) - * t * b * gauss_w5[i] * gauss_w5[j]; - - if (sol!=sol||fabs(sol)==numeric_limits::infinity()){ + 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()) { cout << "->(" << i << "," << j << ")" << endl; - cout << "-|(" << - G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], d, d3) - * t * b * gauss_w5[i] * gauss_w5[j] - << "," << - G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], 0, d3) - * t * b * gauss_w5[i] * gauss_w5[j] - << "," << - G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], d, d3) - * t * b * gauss_w5[i] * gauss_w5[j] - << "," << - G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], 0, d3) - * t * b * gauss_w5[i] * gauss_w5[j] - << ")" << endl; - cout << "||"<< t * b * gauss_w5[i] * gauss_w5[j] <<"||-----------------" << endl; - - cout << "? (" << G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], 0, d3) - << "," << t * b * gauss_w5[i] * gauss_w5[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; @@ -322,9 +365,10 @@ double slpR::calcParIntQY1X1(double b, double d, double t, double v, double d1, } } - if (sol!=sol||sol==numeric_limits::infinity()) - cout << "calcParIntQY1X1: " << sol << " isn't a Number. (" - << b << "," << d << ")(" << t << "," << v << ")(" << d1 << "," << d2 << "," << d3 <<")" << endl; + if (sol != sol || sol == numeric_limits::infinity()) + cout << "calcParIntQY1X1: " << sol << " isn't a Number. (" << b << "," + << d << ")(" << t << "," << v << ")(" << d1 << "," << d2 << "," + << d3 << ")" << endl; return sol; } @@ -333,16 +377,35 @@ double slpR::calcParIntQY1(double b, double d, double t, double v, double d1, double d2, double d3) { //Fall 4 + double sol = 0; + return 0; ///ACHTUNG noch Falsch +} +double slpR::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; + + 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(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) + * gauss_nodes[quad][i].w * gauss_nodes[quad][j].w + * gauss_nodes[quad][k].w * gauss_nodes[quad][l].w; + } + } + } + } - - - return 0; ///ACHTUNG noch Falsch - + return sol; } double slpR::calcOrtIntA(double b, double d, double t, double v, double d1, @@ -356,16 +419,20 @@ double slpR::calcOrtIntQX1X2(double b, double d, double t, double v, double d1, //Fall 2 double sol = 0; - for (int i = 0; i < 5; ++i) { - for (int j = 0; j < 5; ++j) { - sol += G_AY1Y2(-0.5, t + d2, v + d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1) - * b * d * gauss_w5[i] * gauss_w5[j]; - sol -= G_AY1Y2(-0.5, d2, v + d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1) - * b * d * gauss_w5[i] * gauss_w5[j]; - sol -= G_AY1Y2(-0.5, t + d2, d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1) - * b * d * gauss_w5[i] * gauss_w5[j]; - sol += G_AY1Y2(-0.5, d2, d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1) - * b * d * gauss_w5[i] * gauss_w5[j]; + 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; } } @@ -375,11 +442,11 @@ double slpR::calcOrtIntQX1X2(double b, double d, double t, double v, double d1, double slpR::calcParIntO0(double b, double d, double t, double v, double d1, double d2, double d3, double eta) { - return calcParIntA(b,d,t,v,d1,d2,d3); + return calcParIntA(b, d, t, v, d1, d2, d3); } double slpR::calcOrtIntO0(double b, double d, double t, double v, double d1, double d2, double d3, double eta) { - return calcOrtIntA(b,d,t,v,d1,d2,d3); + return calcOrtIntA(b, d, t, v, d1, d2, d3); } double slpR::calcParIntO1(double b, double d, double t, double v, double d1, @@ -446,7 +513,7 @@ double slpR::calcParIntO2(double b, double d, double t, double v, double d1, } // if ((b * b + d * d) * eta < d1 * d1 + d2 * d2 + d3 * d3) { - return calcParIntQX1X2(b, d, t, v, d1, d2, d3); + return calcParIntQX1X2(b, d, t, v, d1, d2, d3); // } else { // cout << "2error"; // return calcParIntA(b, d, t, v, d1, d2, d3); @@ -501,7 +568,7 @@ double slpR::calcParIntO3(double b, double d, double t, double v, double d1, } // if ((min(b, t) * eta <= d1) && (min(d, v) * eta <= d2) ) { - return calcParIntQY1X1(b, d, t, v, d1, d2, d3); + return calcParIntQY1X1(b, d, t, v, d1, d2, d3); // } else { // cout << "3error"; // return calcParIntA(b, d, t, v, d1, d2, d3); @@ -519,36 +586,63 @@ double slpR::calcOrtIntO3(double b, double d, double t, double v, double d1, } } +double slpR::cParO0(double b, double d, double t, double v, double d1, + double d2, double d3, double* eta) { + double tmp = 0; + + // kleine Seite nach vorn + if (b > t) { + tmp = b; + b = t; + t = tmp; + d1 = -d1; + } + + // kleine Seite nach vorn + if (d > v) { + tmp = d; + d = v; + v = tmp; + d2 = -d2; + } + + if ((t * t + v * v) * eta[0] < d1 * d1 + d2 * d2 + d3 * d3) { +// cout << "E"; + return calcParIntQ(b, d, t, v, d1, d2, d3); + } else { + return calcParIntA(b, d, t, v, d1, d2, d3); + } + return 0; +} double slpR::cParO1(double b, double d, double t, double v, double d1, double d2, double d3, double* eta) { - return calcParIntA(b, d, t, v, d1, d2, d3);; + return calcParIntA(b, d, t, v, d1, d2, d3); } - double slpR::cParO2(double b, double d, double t, double v, double d1, double d2, double d3, double* eta) { double tmp = 0; // kleine Seite nach vorn - if(b>t){ + if (b > t) { tmp = b; b = t; t = tmp; - d1=-d1; + d1 = -d1; } // kleine Seite nach vorn - if(d>v){ + if (d > v) { tmp = d; d = v; v = tmp; - d2=-d2; + d2 = -d2; } - if ((b * b + d * d) * eta[0] < d1 * d1 + d2 * d2 + d3 * d3) { + if ((b * b + d * d) < eta[0]* (d1 * d1 + d2 * d2 + d3 * d3)) { // cout << "E"; return calcParIntQX1X2(b, d, t, v, d1, d2, d3); } else { @@ -564,22 +658,22 @@ double slpR::cParO3(double b, double d, double t, double v, double d1, double tmp = 0; // kleine Seite nach vorn - if(b>t){ + if (b > t) { tmp = b; b = t; t = tmp; - d1=-d1; + d1 = -d1; } // kleine Seite nach vorn - if(d>v){ + if (d > v) { tmp = d; d = v; v = tmp; - d2=-d2; + d2 = -d2; } - if ((b * b + d * d) * eta[0] < d1 * d1 + d2 * d2 + d3 * d3) { + if ((b * b + d * d) < eta[0] * (d1 * d1 + d2 * d2 + d3 * d3)) { //kleinere Achse nach vorn if (b * b + t * t > d * d + v * v) { @@ -597,10 +691,10 @@ double slpR::cParO3(double b, double d, double t, double v, double d1, d2 = tmp; } - if (min(b, t) * eta[1] <= d1){ + if (min(b, t) <= eta[1] * d1) { // cout << "A"; return calcParIntQY1X1(b, d, t, v, d1, d2, d3); - }else{ + } else { // cout << "E"; return calcParIntQX1X2(b, d, t, v, d1, d2, d3); } diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp index 6b56825..fae183d 100644 --- a/src/slpRectangle.hpp +++ b/src/slpRectangle.hpp @@ -34,6 +34,7 @@ namespace slpR { 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); @@ -71,6 +72,8 @@ namespace slpR { // sol = cParIntX(b,d,t,v,d1,d2,d3,eta); //Voll Analytisch + double cParO0(double, double, double, double, double, double, double, + double*); double cParO1(double, double, double, double, double, double, double, double*); double cParO2(double, double, double, double, double, double, double,