From ca6c2a19d5658cb76528354c62d40749b3a91fe8 Mon Sep 17 00:00:00 2001 From: treecity Date: Sat, 21 May 2011 16:55:57 +0000 Subject: [PATCH] =?utf8?q?[src]=20Matlab=20g=20G=20und=20F=20entfernt=20[s?= =?utf8?q?rc]=20Art=20der=20Verfeinerung=20Implementiert=20[src]=20refine?= =?utf8?q?=20an=204=20Koo=20pro=20Ele=20angepasst=20[src]=20test=5FsolveE?= =?utf8?q?=20adaptive=20Verfeinerung=20(ob=20verfeinert=20wird=20oder=20ni?= =?utf8?q?cht...=20muss=20noch=20entschieden=20werden)=20[src]=20Beispielg?= =?utf8?q?itter=20in=204Koo=20pro=20Ele=20=C3=BCbersetzt=20mit=20expCOEL?= MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@31 26120e32-c555-405d-b3e1-1f783fb42516 --- src/F_ort.m | 14 ------- src/F_par.m | 22 ----------- src/G00.m | 26 ------------- src/build_A.cpp | 10 ++--- src/exmpl_2DLShape.mat | Bin 555 -> 264 bytes src/exmpl_2DQuad.mat | Bin 243 -> 249 bytes src/exmpl_3DCube.mat | Bin 269 -> 271 bytes src/g0.m | 34 ---------------- src/refineQuad.m | 60 +++++++++++++--------------- src/refineType.m | 25 ++++++++++++ src/test_solveError.m | 86 +++++++++++++++++++++++------------------ 11 files changed, 104 insertions(+), 173 deletions(-) delete mode 100644 src/F_ort.m delete mode 100644 src/F_par.m delete mode 100644 src/G00.m delete mode 100644 src/g0.m create mode 100644 src/refineType.m diff --git a/src/F_ort.m b/src/F_ort.m deleted file mode 100644 index 2c5636b..0000000 --- a/src/F_ort.m +++ /dev/null @@ -1,14 +0,0 @@ -function sol = F_ort(x, y, delta) - - xd = x(1) -delta(1); - xyd=x(2)-y(1)-delta(2); - yd = y(2)+delta(3); - - sol = -G00(0.5,[y(2) x(1)]',[-delta(3) delta(1)],xyd)... - -xd*xyd*G00(-0.5,[x(2) y(2)],[yd -delta(3)],xd)... - +xd*g0(0.5,y(2),-delta(3),norm([xd xyd]))... - -yd*xyd*G00(-0.5,[x(1) x(2)],[delta(1) yd],-yd)... - +yd*g0(0.5,x(1),delta(1),norm([xyd yd])); - -end - diff --git a/src/F_par.m b/src/F_par.m deleted file mode 100644 index c84bb42..0000000 --- a/src/F_par.m +++ /dev/null @@ -1,22 +0,0 @@ -function sol = F_par(x1,x2,y1,y2,d1,d2,d3) -% fprintf('%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f',x1,x2,y1,y2,d1,d2,d3) -% xyd=x-y-delta(1:2); - - sol = (x1-y1-d1)*(x2-y2-d2); - - if(sol~=0) - sol = sol * G00(-0.5,x1,x2,y1+d1,y2+d2,d3); - end - - if((x1-y1-d1)~=0) - sol = sol - (x1-y1-d1)*g0(0.5,x1,y1+d1,sqrt((x2-y2-d2)^2+d3^2)); - end - - if((x2-y2-d2)~=0) - sol = sol -(x2-y2-d2)*g0(0.5,x2,y2+d2,sqrt((x1-y1-d1)^2+d3^2)); - end - - sol = sol + 1/3*((x1-y1-d1)^2+(x2-y2-d2)^2+d3^2)^(1.5); - -end - diff --git a/src/G00.m b/src/G00.m deleted file mode 100644 index b725d43..0000000 --- a/src/G00.m +++ /dev/null @@ -1,26 +0,0 @@ -function sol = G00(p,y1,y2,x1,x2,l) -% fprintf('%.1f | %.1f %.1f | %.1f %.1f | %.1f +',p,x1,x2,y1,y2,l); - sol = 0; - if(p==-1.5) - if(l==0) - sol = -sqrt((y1-x1)^2+(y2-x2)^2)/((y1-x1)*(y2-x2)); - else - sol = sign((y1-x1)*(y2-x2))/(2*abs(l)) ... - * acos(-2*(y1-x1)^2*(y2-x2)^2 ... - /(((y1-x1)^2+l^2)*((y2-x2)^2+l^2))+1); - end - elseif(p==-0.5) - if(l~=0) - sol = 2*p*l^2*G00(p-1,y1,y2,x1,x2,l); - end - if((y1-x1)~=0) - sol = sol + (y1-x1)*g0(p,y2,x2,sqrt((y1-x1)^2+l^2)); - end - if((y2-x2)~=0) - sol = sol + (y2-x2)*g0(p,y1,x1,sqrt((y2-x2)^2+l^2)); - end - sol = sol / (2*p+2); - else - error 'no case for p defined' - end -end diff --git a/src/build_A.cpp b/src/build_A.cpp index 49175fd..4bf70d2 100644 --- a/src/build_A.cpp +++ b/src/build_A.cpp @@ -216,12 +216,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { // for (i = 0; i<3;++i) // d[i] = y0[i]-x0[i]; - printf("(%d,%d)",rx,ry); + //printf("(%d,%d)",rx,ry); if (rx == ry) { - printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb); - printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]); - printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]); - printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]); + //printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb); + //printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]); + //printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]); + //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]); if (rxa == rya) { //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0], // y1[0], y0[1], y3[1], d[0], d[1], d[2]); diff --git a/src/exmpl_2DLShape.mat b/src/exmpl_2DLShape.mat index feb1a9435624f9e922a173ddea746e683980b0e6..655e7bde2958eaedb65001b88e38609ef0253cdb 100644 GIT binary patch delta 31 ncmZ3@(!n&rfydO!(7?*XM8U|w&~Rd)_QV9%i8UM(H?{x(g&zq@ delta 296 zcmeBRTFo-Sfyc}aZCnOj!B$+8ZYq+Fv zfa|I8Awy&5q&P!n;TcB{Ffb^%Fs_@}&{5A1RO|)UumV|wv!Swp6;r04!Z(ITTe2^| z-}liXT|afn>R-7(&8iDcq}1}w3$GrIKKuHO>Vikp0_=R`7&@M+34l!0fSZWo1ZG2J z31cG_#!R&Z9Lmbd%nahprVJM*9+wRViH9%%?V85$=eU3M srHp%bb8iP3tPMAK1F};O845eQYB4jbhunF%!Rf_oK1+56i=`|}09sRM#{d8T diff --git a/src/exmpl_2DQuad.mat b/src/exmpl_2DQuad.mat index eec428172f09e2c930f76b2658f173e25136910e..994a16200fdc7e0fa01c5cd311d3fb83afe68323 100644 GIT binary patch delta 132 zcmV-~0DJ%Q0r>%tH5XxYAWdO;ATlu^F*Q0dH99jjATls9F_BR#kzfUpdzsg;4Lf{}rt;lx1gi3zL|YdG@F7#J8T<~&YL zNZ?>dGE;cga7p0+*HdF>LuG+t=E9=7xSvl@DxO%P;G0;fU}UIZXleySrsfJp28MSdZ>DqfnW;8;+kU}&IVXli9@X=P}kU}RuuI5ALrVgl>L8jgB<1_p+TIggVQ z5;_=?%oLtATv9l|^>j0{p|OB5lcBkQ`t#@FOH-byPFo_@eASR~qXu&XqjH3xva*Ze z1ja=Se|1= sqrt(Ct(2,:).^2+Ct(4,:).^2)); +% REF(t3) = 3; +% t4 = (tau*abs(Ct(4,:)) >= sqrt(Ct(2,:).^2+Ct(3,:).^2)); +% REF(t4) = 4; +REF(~(t4+t3+t1)) = 2; +end + diff --git a/src/test_solveError.m b/src/test_solveError.m index 4eca864..9730b85 100644 --- a/src/test_solveError.m +++ b/src/test_solveError.m @@ -1,51 +1,62 @@ +%% Laden des Gitters +load exmpl_2DLShape %Energienorm sollte gegen 8.28466 konvergieren -load exmpl_2DLShape - +%% Rotationen einzelnen von Elementen (zum Testen) % coordinates = coordinates(:,[ 1 3 2]); -elements(2,:) = elements(2,[ 3 4 1 2]); -elements(3,:) = elements(3,[ 3 2 1 4]); +% elements(2,:) = elements(2,[ 3 4 1 2]); +% elements(3,:) = elements(3,[ 3 2 1 4]); +%% Verfeinerungen (isotrop) +for i = 1:0 + [coordinates,elements]=refineQuad(coordinates,elements,2*ones(1,size(elements,1))); +end % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1))); % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1))); % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1))); % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1))); -% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1))); -% build_A(coordinates,elements); -% A = build_A2(coordinates,elements); -% sum(sum (A-A')); -% b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2)); -% x = A\b; -% xe = x'*A*x +for ref = 1:4 -A = build_A(coordinates,elements) -sum(sum (A-A')); -b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2)); -x = A\b; -xe = x'*A*x + %% Testen mit Quadratur... und mex_G00 + % A = build_A2(coordinates,elements); + % sum(sum (A-A')); + % b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2)); + % x = A\b; + % xe = x'*A*x -% [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1))); -% A_fine = build_A2(coordinates_fine,elements_fine) -% sum(sum(A_fine-A_fine')) -% b = sqrt(sum(quad_norm(coordinates_fine,elements_fine,'w').^2,2)); -% x_fine = A_fine\b; -% xe_fine = x_fine'*A_fine*x_fine; -% -% x_interpol = zeros(size(elements_fine,1),1); -% for i = 1:size(elements,1) -% x_interpol(f2s(i,:)) = x(i); -% end -% xe_interpol = x_interpol'*A_fine*x_interpol; -% -% [x_fine x_interpol x_fine-x_interpol] -% -% % l = x_fine(f2s(1,:)); -% % T4 = [1 1 1 1; 1 -1 1 -1; 1 1 -1 -1;1 -1 -1 1]; -% % l -% % T4*l*1/4 -% -% [xe xe_fine xe_interpol] + %% Testen mit MEX -> build_A + A = build_A(coordinates,elements); + sum(sum (A-A')); + b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2)); + x = A\b; + xe = x'*A*x; + + %% isotrope Verfeinerung und Aufbauen + [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1))); + A_fine = build_A(coordinates_fine,elements_fine); + sum(sum(A_fine-A_fine')); + b = sqrt(sum(quad_norm(coordinates_fine,elements_fine,'w').^2,2)); + x_fine = A_fine\b; + xe_fine = x_fine'*A_fine*x_fine; + + + x_interpol = zeros(size(elements_fine,1),1); + for i = 1:size(elements,1) + x_interpol(f2s(i,:)) = x(i); + end + xe_interpol = x_interpol'*A_fine*x_interpol; + + % [x_fine x_interpol x_fine-x_interpol] + + % Wie soll verfeinert werden + marked = refineType(x_fine(f2s)',0.3); + + [xe xe_fine xe_interpol] + + [coordinates, elements] = refineQuad(coordinates,elements,marked); + +end % re_times = 6; % r_type = 1; @@ -68,4 +79,3 @@ xe = x'*A*x % legend('x^TAx') % xlabel('Elemente') % ylabel('xNorm') - -- 2.47.3