[src] 4fach Gauss durch doppelgauss mit G00 ersetzt

[src] verzweifelte Fehlersuche... Matrix stimmt mit Gauss Matrix überein...

git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@26 26120e32-c555-405d-b3e1-1f783fb42516

diff --git a/src/build_A.cpp b/src/build_A.cpp
index 9281525d30711ff6293c9b96d04f9341fef744ea..60504802cdeda10b9594b8a88db745a3f436c16b 100644 (file)
--- 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[]) {
 \r
                        if (rx = ry) {\r
                                if (rx = 2) {\r
-                                       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],\r
-                                                       y2[0], y1[1], y3[1], d[0], d[1], d[2]);\r
+                                       //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],\r
+                                       //              y2[0], y1[1], y3[1], d[0], d[1], d[2]);\r
                                        //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
                                        tmp = quadInt(F_par, x1[0], x2[0], x1[1], x3[1], y1[0],\r
                                                        y2[0], y1[1], y3[1], d[0], d[1], d[2]);\r

diff --git a/src/build_A2.m b/src/build_A2.m
index 9ca4079ff320f608302e298bb870470d137e62dd..4956416d102d37a5e2109d8d2389cda33168913c 100644 (file)
--- 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 17fed8f3f84ddd6643506bb643109eba8f14d945..757d135497a9dccf009e092b136f651c3225bccb 100644 (file)
--- 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 0049235eeb3129ab8260c5c754ed8a177f35c3b7..366ea259b910d9e389579c4a8a2c5b953759342b 100644 (file)
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -1,40 +1,46 @@
 \r
 \r
 load exmpl_2DLShape\r
-\r
-\r
+%[coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
 % [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
-\r
+% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
+ [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
 % build_A(coordinates,elements);\r
 \r
-A = build_A2(coordinates,elements)\r
-sum(sum (A-A'))\r
+A = build_A2(coordinates,elements);\r
+sum(sum (A-A'));\r
+b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2))\r
+x = A\b;\r
+xe = x'*A*x\r
+\r
+A = build_A(coordinates,elements);% .* (ones(size(elements,1))+eye(size(elements,1)));\r
+sum(sum (A-A'));\r
 b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));\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_A2(coordinates_fine,elements_fine)\r
-sum(sum(A_fine-A_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);\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 xe_interpol]\r
+% [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
+% A_fine = build_A2(coordinates_fine,elements_fine)\r
+% sum(sum(A_fine-A_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);\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 xe_interpol]\r
 \r
 % re_times = 6;\r
 % r_type = 1;\r