[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

diff --git a/src/F_par.m b/src/F_par.m
index b3e5ebb289ca70cf584af16a7e4997c2675290be..c84bb429e6f1d7cd9844a9e5fd8c45ae71009d08 100644 (file)
--- 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 (file)
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 15714196bb4fae7a5d8a50d4917d297a22dfbbaa..023fa2d97452df0fd8e42f16547a8ab6107042ac 100644 (file)
--- 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 848a7069b54b63e8c0042a16686be0d6bb108eed..e0c2551038afb4395232038b4b8fadc8772e29c2 100644 (file)
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -2,23 +2,31 @@
 \r
 load exmpl_2DLShape\r
 \r
-A = build_A(coordinates,elements);\r
+A = build_A2(coordinates,elements);\r
 b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));\r
-x = A\b;\r
+x = A\b\r
 xe = x'*A*x;\r
 \r
 [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
-A_fine = build_A(coordinates_fine,elements_fine);\r
+A_fine = build_A2(coordinates_fine,elements_fine);\r
 b = sqrt(sum(quad_norm(coordinates_fine,elements_fine,'w').^2,2));\r
 x_fine = A_fine\b;\r
 xe_fine = x_fine'*A_fine*x_fine;\r
 \r
+x_interpol = zeros(size(elements_fine,1),1);\r
+for i = 1:size(elements,1)\r
+    x_interpol(f2s(i,:)) = x(i)/4;\r
+end\r
+xe_interpol = x_interpol'*A_fine*x_interpol;\r
+\r
+[x_fine x_interpol x_fine-x_interpol]\r
+\r
 l = x_fine(f2s(1,:));\r
 T4 = [1 1 1 1; 1 -1 1 -1; 1 1 -1 -1;1 -1 -1 1];\r
 l\r
 T4*l*1/4\r
 \r
-[xe xe_fine]\r
+[xe xe_fine xe_interpol]\r
 \r
 % re_times = 6;\r
 % r_type = 1;\r