From eeebc09be1d32c317b1f16865357cbd60d20cdf8 Mon Sep 17 00:00:00 2001 From: treecity Date: Mon, 16 May 2011 15:45:52 +0000 Subject: [PATCH] =?utf8?q?[src]=204fach=20Gauss=20durch=20doppelgauss=20mi?= =?utf8?q?t=20G00=20ersetzt=20[src]=20verzweifelte=20Fehlersuche...=20Matr?= =?utf8?q?ix=20stimmt=20mit=20Gauss=20Matrix=20=C3=BCberein...?= 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@26 26120e32-c555-405d-b3e1-1f783fb42516 --- src/build_A.cpp | 4 +-- src/build_A2.m | 12 ++++----- src/surfDoubleQuad.m | 50 +++++++++++++++++++++++++++++++------ src/test_solveError.m | 58 ++++++++++++++++++++++++------------------- 4 files changed, 83 insertions(+), 41 deletions(-) diff --git a/src/build_A.cpp b/src/build_A.cpp index 9281525..6050480 100644 --- a/src/build_A.cpp +++ b/src/build_A.cpp @@ -114,8 +114,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (rx = ry) { if (rx = 2) { - printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x1[0], x2[0], x1[1], x3[1], y1[0], - y2[0], y1[1], y3[1], d[0], d[1], d[2]); + //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x1[0], x2[0], x1[1], x3[1], y1[0], + // y2[0], y1[1], y3[1], d[0], d[1], d[2]); //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]); tmp = quadInt(F_par, x1[0], x2[0], x1[1], x3[1], y1[0], y2[0], y1[1], y3[1], d[0], d[1], d[2]); diff --git a/src/build_A2.m b/src/build_A2.m index 9ca4079..4956416 100644 --- a/src/build_A2.m +++ b/src/build_A2.m @@ -4,11 +4,11 @@ function A = build_A2(coordinates,elements) A1 = zeros(N); % untere schranke s t obere schranke k l - intF = @(f,a,b,c,d,r,t,u,v)... - f(b,d,t,v)-f(b,d,t,u)-f(b,d,r,v)+f(b,d,r,u)... - -f(b,c,t,v)+f(b,c,t,u)+f(b,c,r,v)-f(b,c,r,u)... - -f(a,d,t,v)+f(a,d,t,u)+f(a,d,r,v)-f(a,d,r,u)... - +f(a,c,t,v)-f(a,c,t,u)-f(a,c,r,v)+f(a,c,r,u); +% intF = @(f,a,b,c,d,r,t,u,v)... +% f(b,d,t,v)-f(b,d,t,u)-f(b,d,r,v)+f(b,d,r,u)... +% -f(b,c,t,v)+f(b,c,t,u)+f(b,c,r,v)-f(b,c,r,u)... +% -f(a,d,t,v)+f(a,d,t,u)+f(a,d,r,v)-f(a,d,r,u)... +% +f(a,c,t,v)-f(a,c,t,u)-f(a,c,r,v)+f(a,c,r,u); tic for j = 1:N for k = 1:N @@ -16,7 +16,7 @@ function A = build_A2(coordinates,elements) ek = coordinates(elements(k,:)',:); d = ek(1,:) - ej(1,:); % d = ones(1,3); - A1(j,k) = surfDoubleQuad(@(x1,x2,y1,y2)1/sqrt((x1-y1-d(1)).^2+(x2-y2-d(2)).^2+d(3).^2)... + A1(j,k) = surfDoubleQuad(@(x1,x2,y1,y2) 1/sqrt((x1-y1-d(1)).^2+(x2-y2-d(2)).^2+d(3).^2)... ,ej(1,1),ej(2,1),ej(1,2), ej(3,2),ek(1,1), ek(2,1),ek(1,2), ek(3,2),4); % A1(j,k) = intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))... diff --git a/src/surfDoubleQuad.m b/src/surfDoubleQuad.m index 17fed8f..757d135 100644 --- a/src/surfDoubleQuad.m +++ b/src/surfDoubleQuad.m @@ -45,24 +45,60 @@ % end +% function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w) +% +% [x1n x1w] = gauss(w,a,b); +% [x2n x2w] = gauss(w,c,d); +% [y1n y1w] = gauss(w,r,t); +% [y2n y2w] = gauss(w,u,v); +% +% in = 0; +% for i=1:length(y2n) +% for j=1:length(y1n) +% for k=1:length(x2n) +% for l = 1:length(x1n) +% +% in = in + y2w(i) * y1w(j) * x2w(k) * x1w(l) * f(x1n(l),x2n(k),y1n(j),y2n(i)); +% +% end +% end +% end +% end +% +% end + function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w) [x1n x1w] = gauss(w,a,b); [x2n x2w] = gauss(w,c,d); - [y1n y1w] = gauss(w,r,t); - [y2n y2w] = gauss(w,u,v); + + + d1=r-a; + d2=u-c; +% sum(x1w)a + x1n=x1n-d1; x2n=x2n-d2; +% c +% d + + p = -0.5; in = 0; - for i=1:length(y2n) - for j=1:length(y1n) + for k=1:length(x2n) for l = 1:length(x1n) - in = in + y2w(i) * y1w(j) * x2w(k) * x1w(l) * f(x1n(l),x2n(k),y1n(j),y2n(i)); + [my(1) mk(1)] = mex_G00(p,r,u,x1n(l),x2n(k),0); + [my(2) mk(2)] = mex_G00(p,t,v,x1n(l),x2n(k),0); + [my(3) mk(3)] = mex_G00(p,t,u,x1n(l),x2n(k),0); + [my(4) mk(4)] = mex_G00(p,r,v,x1n(l),x2n(k),0); + + +% my_sol = my(1)+my(2)-my(3)-my(4); + + in = in + x2w(k) * x1w(l) * (my(1)+my(2)-my(3)-my(4)); end end - end - end + in end \ No newline at end of file diff --git a/src/test_solveError.m b/src/test_solveError.m index 0049235..366ea25 100644 --- a/src/test_solveError.m +++ b/src/test_solveError.m @@ -1,40 +1,46 @@ load exmpl_2DLShape - - +%[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')) +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 + +A = build_A(coordinates,elements);% .* (ones(size(elements,1))+eye(size(elements,1))); +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] +% [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] % re_times = 6; % r_type = 1; -- 2.47.3