From afb01996bf6efd3fb3386c7b302b0bbe9212711b Mon Sep 17 00:00:00 2001 From: treecity Date: Tue, 8 Nov 2011 12:49:57 +0000 Subject: [PATCH] [src] weitere Anpassungen an die test_calc Files [src] erste Versuche den QuadAchse NAN Fehler zu finden git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@66 26120e32-c555-405d-b3e1-1f783fb42516 --- src/slpRectangle.cpp | 43 ++++++++++++++++++++++++++---- src/test_calcInt1.m | 2 +- src/test_calcInt2.m | 60 ++++++++++++++++++++++++++++++++--------- src/test_calcInt3.m | 63 ++++++++++++++++++++++++++++++++++---------- 4 files changed, 135 insertions(+), 33 deletions(-) diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp index 3dbe4c5..84d8175 100644 --- a/src/slpRectangle.cpp +++ b/src/slpRectangle.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include //#include @@ -97,15 +98,20 @@ double G_QY1Y2(double p, double y1, double y2, double x1, double x2, double l) { double G_AY2X2(double y1, double y2, double x1, double x2, double d3) { - // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3); +// cout << "(" << y1 << "," << y2 << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl; double hlp = ((x1 - y1) * (x1 - y1) + (x2 - y2) * (x2 - y2) + d3 * d3); double sol = sqrt(hlp); + if ((x2 - y2) != 0) sol += (x2 - y2) * log(sqrt(hlp) - (x2 - y2)); + if (sol!=sol||sol==numeric_limits::infinity()) + cerr << "G_AY2X2: " << sol << " isn't a Number. (" + << y1 << "," << y2 << ")(" << x1 << "," << x2 << ")(" << d3 << ")" << endl; + return sol; } @@ -274,19 +280,46 @@ double calcParIntQY1X1(double b, double d, double t, double v, double d1, 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]; + += 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]; + -= 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||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 << endl; + return sol; + } + } } + if (sol!=sol||sol==numeric_limits::infinity()) + cout << "calcParIntQY1X1: " << sol << " isn't a Number. (" + << b << "," << d << ")(" << t << "," << v << ")(" << d1 << "," << d2 << "," << d3 <<")" << endl; + return sol; } diff --git a/src/test_calcInt1.m b/src/test_calcInt1.m index ad6da47..b342959 100644 --- a/src/test_calcInt1.m +++ b/src/test_calcInt1.m @@ -1 +1 @@ - coo = @(h,diff)[0 0 0;1 0 0; 1 1 0;0 1 0; 0 0 0 ; 0+h 0 0; 0+h h 0 ; 0 h 0]+[zeros(4,3);repmat(diff,4,1)]; elements=[1 2 3 4;5 6 7 8]; neigh = zeros(2,8); dat = []; diff = [2 0 0]; %% Laage Übersicht figure(1) h = 1; coordinates=coo(h,diff) current = coordinates(elements(1,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'g'); hold on current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'b'); h = h/2; coordinates=coo(h,diff); current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'y'); hold off legend('Element1(h)','Element2(h)','Element2(h/2)'); title('Laage der Elemente'); %% Integrale bei kleiner werdenden Elementen h = 2; for I = 1:50 h = h/2; coordinates=coo(h,diff); A0 = mex_build_AU(coordinates,elements,0,0); A2 = mex_build_AU(coordinates,elements,1,2); A1 = mex_build_AU(coordinates,elements,1,1); I dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)]; dat(dat<0)=0; end figure(2) loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4))); legend('Analytisch','Quad Element','Element vertauschen','location','southeast'); xlabel('Elementgroesse (kürzeste Seite)'); ylabel('Integral'); title('Integral bei kleiner werdenden Element'); %% Netzverfeinerung ab einem Wert h = (1/2)^25; coordinates=coo(h,diff); A_loadMesh(coordinates,elements,neigh); datA=[]; for I = 1:5 datA(I,:) = A_stepIso(5,[0 1 3 2]); % datA(I,:) = A_stepAniso(6,[0 1 3 2],1,0); I end figure(3) loglog(datA(2:end,1),datA(2:end,3),datA(2:end,1),datA(2:end,6),datA(2:end,1),datA(2:end,12),datA(2:end,1),datA(2:end,9)) legend('Analytisch','Element vertauschen','Quad Element','Quad Achse','location','northeast'); xlabel('Anzahl der Elemente'); ylabel('mu Schaetzer'); title('mu Schaetzer mit "schlechtem" Startnetz'); figure(4) loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,7),datA(:,1),datA(:,13),datA(:,1),datA(:,10)) legend('Analytisch','Element vertauschen','Quad Element','Quad Achse','location','southeast'); xlabel('Anzahl der Elemente'); ylabel('EnergieNorm ^2 '); title('EnergieNorm ^2 mit "schlechtem" Startnetz'); datA \ No newline at end of file + coo = @(h,diff)[0 0 0;1 0 0; 1 1 0;0 1 0; 0 0 0 ; 0+h 0 0; 0+h h 0 ; 0 h 0]+[zeros(4,3);repmat(diff,4,1)]; elements=[1 2 3 4;5 6 7 8]; neigh = zeros(2,8); dat = []; diff = [2 0 0]; %% Laage Übersicht figure(1) h = 1; coordinates=coo(h,diff) current = coordinates(elements(1,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'g'); hold on current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'b'); h = h/2; coordinates=coo(h,diff); current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'y'); hold off legend('Element1(h)','Element2(h)','Element2(h/2)'); title('Laage der Elemente'); %% Integrale bei kleiner werdenden Elementen h = 2; for I = 1:50 h = h/2; coordinates=coo(h,diff); A0 = mex_build_AU(coordinates,elements,0,0); A2 = mex_build_AU(coordinates,elements,1,2); A1 = mex_build_AU(coordinates,elements,1,1); I dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)]; dat(dat<0)=0; end figure(2) loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4))); legend('Analytisch','Quad Element','Element vertauschen','location','best'); xlabel('Elementgroesse (kürzeste Seite)'); ylabel('Integral'); title('Integral bei kleiner werdenden Element'); %% Netzverfeinerung mit schlechtem Netz h = (1/2)^25; coordinates=coo(h,diff); A_loadMesh(coordinates,elements,neigh); datA=[]; for I = 1:5 datA(I,:) = A_stepIso(1,[0 3 2]); % datA(I,:) = A_stepAniso(1,[0 3 2],1,0); I end figure(3) loglog(datA(:,1),datA(:,3),datA(:,1),datA(:,9),datA(:,1),datA(:,6)) legend('Analytisch','Quad Element','Quad Achse','location','best'); xlabel('Anzahl der Elemente'); ylabel('mu Schaetzer'); title('mu Schaetzer mit "schlechtem" Startnetz'); figure(4) loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7)) legend('Analytisch','Quad Element','Quad Achse','location','best'); xlabel('Anzahl der Elemente'); ylabel('EnergieNorm ^2 '); title('EnergieNorm ^2 mit "schlechtem" Startnetz'); datA \ No newline at end of file diff --git a/src/test_calcInt2.m b/src/test_calcInt2.m index c86ccfc..34a9e2d 100644 --- a/src/test_calcInt2.m +++ b/src/test_calcInt2.m @@ -1,15 +1,16 @@ +coo = @(h,diff)[0 0 0;1 0 0; 1 1 0;0 1 0; 0 0 0 ; 0+h 0 0; 0+h 1 0 ; 0 1 0]+[zeros(4,3);repmat(diff,4,1)]; elements=[1 2 3 4;5 6 7 8]; neigh = zeros(2,8); - dat = []; -diff = 10; +diff = [5 0 0]; -h = 1/2; -coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; diff 0 0; diff+h 0 0; diff+h 1 0 ; diff 1 0]; +%% Laage Übersicht figure(1) +h = 1; +coordinates=coo(h,diff) current = coordinates(elements(1,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'g'); hold on @@ -17,30 +18,63 @@ current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'b'); h = h/2; -coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; diff 0 0; diff+h 0 0; diff+h 1 0 ; diff 1 0]; +coordinates=coo(h,diff); current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'y'); hold off -legend('Element1(h/2)','Element2(h/2)','Element2(h/4)') +legend('Element1(h)','Element2(h)','Element2(h/2)'); +title('Laage der Elemente'); + +%% Integrale bei kleiner werdenden Elementen h = 2; for I = 1:50 h = h/2; -coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; diff 0 0; diff+h 0 0; diff+h 1 0 ; diff 1 0]; + coordinates=coo(h,diff); A0 = mex_build_AU(coordinates,elements,0,0); A2 = mex_build_AU(coordinates,elements,1,2); - A3 = mex_build_AU(coordinates,elements,1,3); + A1 = mex_build_AU(coordinates,elements,1,3); I - dat(I,1:4) = [h A0(1,2) A2(1,2) A3(1,2)]; -% dat(dat<0)=0; + dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)]; + dat(dat<0)=0; end figure(2) loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4))); -legend('Analytisch','Quad Element','Quad Achse'); +legend('Analytisch','Quad Element','Quad Achse','location','best'); +xlabel('Elementgroesse (kürzeste Seite)'); +ylabel('Integral'); +title('Integral bei kleiner werdenden Element'); + +%% Netzverfeinerung mit schlechtem Netz + +h = (1/2)^25; +coordinates=coo(h,diff); + +A_loadMesh(coordinates,elements,neigh); +datA=[]; +for I = 1:5 + datA(I,:) = A_stepIso(1,[0 3 2]); +% datA(I,:) = A_stepAniso(1,[0 1 3 2],1,0); + I +end +figure(3) +loglog(datA(:,1),datA(:,3),datA(:,1),datA(:,9),datA(:,1),datA(:,6)) + +legend('Analytisch','Quad Element','Quad Achse','location','best'); +xlabel('Anzahl der Elemente'); +ylabel('mu Schaetzer'); +title('mu Schaetzer mit "schlechtem" Startnetz'); + +figure(4) +loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7)) + +legend('Analytisch','Quad Element','Quad Achse','location','best'); +xlabel('Anzahl der Elemente'); +ylabel('EnergieNorm ^2 '); +title('EnergieNorm ^2 mit "schlechtem" Startnetz'); -xlabel('Elementgroesse'); -ylabel('Integral'); \ No newline at end of file +datA diff --git a/src/test_calcInt3.m b/src/test_calcInt3.m index f83514e..d530636 100644 --- a/src/test_calcInt3.m +++ b/src/test_calcInt3.m @@ -1,13 +1,16 @@ +coo = @(h,diff)[0 0 0;h 0 0; h 1 0;0 1 0; 0 0 0 ; 0+h 0 0; 0+h 1 0 ; 0 1 0]+[zeros(4,3);repmat(diff,4,1)]; elements=[1 2 3 4;5 6 7 8]; neigh = zeros(2,8); - dat = []; -h = 1/2; -coordinates=[0 0 0;h 0 0; h 1 0;0 1 0; 2 0 0; 2+h 0 0; 2+h 1 0 ; 2 1 0]; +diff = [5 0 0]; + +%% Laage Übersicht figure(1) +h = 1; +coordinates=coo(h,diff) current = coordinates(elements(1,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'g'); hold on @@ -15,33 +18,65 @@ current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'b'); h = h/2; -coordinates=[0 0 0;h 0 0; h 1 0;0 1 0; 2 0 0; 2+h 0 0; 2+h 1 0 ; 2 1 0]; +coordinates=coo(h,diff); current = coordinates(elements(1,[1:4,1])',:); -fill3(current(:,1),current(:,2),current(:,3),'y'); -hold on -current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'r'); +current = coordinates(elements(2,[1:4,1])',:); +fill3(current(:,1),current(:,2),current(:,3),'y'); hold off -legend('Element1(h/2)','Element2(h/2)','Element1(h/4)','Element2(h/4)') +legend('Element1(h)','Element2(h)','Element1(h/2)','Element2(h/2)'); +title('Laage der Elemente'); + +%% Integrale bei kleiner werdenden Elementen h = 2; for I = 1:50 h = h/2; -coordinates=[0 0 0;h 0 0; h 1 0;0 1 0; 2 0 0; 2+h 0 0; 2+h 1 0 ; 2 1 0]; + coordinates=coo(h,diff); A0 = mex_build_AU(coordinates,elements,0,0); A2 = mex_build_AU(coordinates,elements,1,2); - A3 = mex_build_AU(coordinates,elements,1,3); + A1 = mex_build_AU(coordinates,elements,1,3); I - dat(I,1:4) = [h A0(1,2) A2(1,2) A3(1,2)]; + dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)]; dat(dat<0)=0; end figure(2) loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4))); -legend('Analytisch','Quad Element','Quad Achse'); +legend('Analytisch','Quad Element','Quad Achse','location','best'); +xlabel('Elementgroesse (kürzeste Seite)'); +ylabel('Integral'); +title('Integral bei kleiner werdenden Element'); + +%% Netzverfeinerung mit schlechtem Netz + +h = (1/2)^25; +coordinates=coo(h,diff); + +A_loadMesh(coordinates,elements,neigh); +datA=[]; +for I = 1:5 + datA(I,:) = A_stepIso(1,[0 3 2]); +% datA(I,:) = A_stepAniso(1,[0 3 2],1,0); + I +end +figure(3) +loglog(datA(:,1),datA(:,3),datA(:,1),datA(:,9),datA(:,1),datA(:,6)) + +legend('Analytisch','Quad Element','Quad Achse','location','best'); +xlabel('Anzahl der Elemente'); +ylabel('mu Schaetzer'); +title('mu Schaetzer mit "schlechtem" Startnetz'); + +figure(4) +loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7)) + +legend('Analytisch','Quad Element','Quad Achse','location','best'); +xlabel('Anzahl der Elemente'); +ylabel('EnergieNorm ^2 '); +title('EnergieNorm ^2 mit "schlechtem" Startnetz'); -xlabel('Elementgroesse'); -ylabel('Integral'); \ No newline at end of file +datA -- 2.47.3