From 988f026e9b274275e6418454aa1ce1eb1fce10e0 Mon Sep 17 00:00:00 2001 From: treecity Date: Mon, 2 May 2011 17:28:35 +0000 Subject: [PATCH] [src] auch wieder auf Matlab build umgestiegen um den Fehler schneller zu finden git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@18 26120e32-c555-405d-b3e1-1f783fb42516 --- src/F_par.m | 11 ----------- src/build_A2.m | 26 ++++++++++++++++++++++++++ src/g0.m | 2 +- src/test_solveError.m | 16 ++++++++++++---- 4 files changed, 39 insertions(+), 16 deletions(-) create mode 100644 src/build_A2.m diff --git a/src/F_par.m b/src/F_par.m index b3e5ebb..c84bb42 100644 --- a/src/F_par.m +++ b/src/F_par.m @@ -18,16 +18,5 @@ function sol = F_par(x1,x2,y1,y2,d1,d2,d3) sol = sol + 1/3*((x1-y1-d1)^2+(x2-y2-d2)^2+d3^2)^(1.5); - -% sol = (x1-y1-d1)*(x2-y2-d2)*G00(-0.5,x1,x2,y1+d1,y2+d2,d3)... -% -(x1-y1-d1)*g0(0.5,x1,y1+d1,sqrt((x2-y2-d2)^2+d3^2))... -% -(x2-y2-d2)*g0(0.5,x2,y2+d2,sqrt((x1-y1-d1)^2+d3^2))... -% +1/3*((x1-y1-d1)^2+(x2-y2-d2)^2+d3^2)^(1.5); - -% sol = xyd(1)*xyd(2)*G00(-0.5,x,y+delta(1:2),delta(3))... -% -xyd(1)*g0(0.5,x(1),y(1)+delta(1),norm([xyd(2) delta(3)]))... -% -xyd(2)*g0(0.5,x(2),y(2)+delta(2),norm([xyd(1) delta(3)]))... -% +1/3*(xyd(1)^2+xyd(2)^2+delta(3)^2)^(1.5); - end diff --git a/src/build_A2.m b/src/build_A2.m new file mode 100644 index 0000000..f6765a9 --- /dev/null +++ b/src/build_A2.m @@ -0,0 +1,26 @@ +function A = build_A2(coordinates,elements) + N = size(elements,1); + + A1 = zeros(N); + + % untere schranke s t obere schranke k l + intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ... + f(k1,k2,l1,l2)-f(k1,k2,l1,t2)-f(k1,k2,t1,l2)+f(k1,k2,t1,t2)... + -f(k1,s2,l1,l2)+f(k1,s2,l1,t2)+f(k1,s2,t1,l2)-f(k1,s2,t1,t2)... + -f(s1,k2,l1,l2)+f(s1,k2,l1,t2)+f(s1,k2,t1,l2)-f(s1,l2,t1,t2)... + +f(s1,s2,l1,l2)-f(s1,0,l1,t2)-f(s1,s2,t1,l2)+f(s1,s2,t1,t2); + tic + for j = 1:N + for k = 1:N + ej = coordinates(elements(j,:)',:); + ek = coordinates(elements(k,:)',:); + d = ek(1,:) - ej(1,:); +% d = ones(1,3); + A1(j,k) = intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))... + ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2)); + end + end + + A = A1/(4*pi); + +end \ No newline at end of file diff --git a/src/g0.m b/src/g0.m index 1571419..023fa2d 100644 --- a/src/g0.m +++ b/src/g0.m @@ -16,7 +16,7 @@ function sol = g0(p, y, x, l) elseif(p==-1.5) sol = yx/(l^2)*((yx^2 + l^2)^-0.5); else - sol = yx*((yx^2 + l^2)^p)+2*p*l^2*g0(p-1,y,x,l); + sol = yx*((yx^2 + l^2)^p)+2*p*l^2*g0(p-1,y,x,l)/(2*p+1); end else if(p==-0.5) diff --git a/src/test_solveError.m b/src/test_solveError.m index 848a706..e0c2551 100644 --- a/src/test_solveError.m +++ b/src/test_solveError.m @@ -2,23 +2,31 @@ load exmpl_2DLShape -A = build_A(coordinates,elements); +A = build_A2(coordinates,elements); b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2)); -x = A\b; +x = A\b xe = x'*A*x; [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1))); -A_fine = build_A(coordinates_fine,elements_fine); +A_fine = build_A2(coordinates_fine,elements_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)/4; +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 xe_fine xe_interpol] % re_times = 6; % r_type = 1; -- 2.47.3