[src] Quadratur wieder aktiviert

[src] Quadratur wieder angepasst

diff --git a/src/A_plots.m b/src/A_plots.m
index 0b603b726371a5312df282a6512933a6be4ff125..a612e6caca26955a43c97121d0d76f4f46cefebf 100644 (file)
--- a/src/A_plots.m
+++ b/src/A_plots.m
@@ -55,7 +55,7 @@ for i = 1:length(files)
           ['min hmax/max hmax ' l0 l1{i}]...
           ['min hmin/hmax ' l0 l1{i}]...
           }';
-      sym = {sym{:} type2sym{data(1,[2+(i-1)*rows])}}'
+      sym = {sym{:} type2sym{data(1,[2+(i-1)*rows])}}';
    end
    
 end
@@ -67,7 +67,7 @@ if step<1
     disp ('Error: No Data to show.')
 else
 
-    data((end-8):end,[1 [3 4 5]])
+%     data((end-8):end,[1 [3 4 5]])
     
 %     sol = interp1(1./X((round(1)):(end),4)',G_D((round(1)):(end),4)',0,'spline')
 %     sol = 8.28466; % LShape

diff --git a/src/A_step.m b/src/A_step.m
index 71657c138f150460a57d2bcfa09bcef976e1c65d..038d47f4de7927583e62bfdf3bc532169c4fef7f 100644 (file)
--- a/src/A_step.m
+++ b/src/A_step.m
@@ -44,8 +44,8 @@ time = zeros(1,3);
   %data -> ErgebnisMatrix
   data = size(G_E,1);
   %save_* -> ZwischenSpeicherung
-%  save_A = {};
-%  save_x = {};
+ save_A = {};
+ save_x = {};
   
   %Alle MatrixBrechenungsArten mit dem selben Netz berechnen
   for i = 1:length(type)
@@ -90,13 +90,13 @@ time = zeros(1,3);
     xe_fine = x_fine'*A_fine*x_fine;
     xe = x'*A*x;
     
-%    save_A{i} = A_fine;
-%    save_x{i} = x_fine;
+   save_A{i} = A_fine;
+   save_x{i} = x_fine;
     
     
     
     data = [data type(i) sqrt(sum(tmu)) sqrt(eta) xe_fine sqrt(sum(mu))...
-        min(hmin)/max(hmax) min(hmax)/max(hmax) sqrt(sum((hmin./hmax).^2))];
+        min(hmin)/max(hmax) min(hmax)/max(hmax) min(hmin./hmax)];
   end
   time(2) = toc;
   
@@ -141,8 +141,8 @@ time = zeros(1,3);
      out = varargin{1};
    end
   typeN = int2str(type);
-%  save (['meshSave/fine_' out typeN(typeN~=' ') '_' int2str(size(G_T,1))],...
-%      'coo_fine', 'ele_fine','neigh_fine','f2s','time','data','save_A','save_x')
+ save (['meshSave/fine_' out typeN(typeN~=' ') '_' int2str(size(G_T,1))],...
+     'coo_fine', 'ele_fine','neigh_fine','f2s','time','data','save_A','save_x')
 
 %   clear 'coo_fine' 'ele_fine' 'neigh_fine' 'f2s' 
 %   plotShape(G_C,G_E,'');

diff --git a/src/build_A2.m b/src/build_A2.m
new file mode 100644 (file)
index 0000000..c28c949
--- /dev/null
+++ b/src/build_A2.m
@@ -0,0 +1,44 @@
+function A = build_A2(coordinates,elements)
+    N = size(elements,1);
+
+     A1 = zeros(N);
+%     A1 = mex_build_V(coordinates,elements,1,1);
+
+    % 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);
+    
+%     matlabpool
+    
+    for j = 1:N
+        for k = j+1:N
+            ej = coordinates(elements(j,[1,2,4])',:);
+            ek = coordinates(elements(k,[1,2,4])',:);
+            
+%             d = ek(1,:) - ej(1,:)
+            
+%             ej = ej - repmat(ej(1,:),3,1);
+%             ek = ek - repmat(ek(1,:),3,1);
+
+            d = zeros(1,3);
+%             if(j~=k)
+            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(k,j) = A1(j,k);
+%             else
+                
+%             end
+
+%             A1(j,k) = intF(@(x1,x2,y1,y2) mex_Fpar(x1,x2,y1,y2,d(1),d(2),d(3))...
+%                 ,ej(1,1),ej(2,1),ej(1,2), ej(3,2),ek(1,1), ek(2,1),ek(1,2), ek(3,2));
+        end
+    end
+    
+%     matlabpool close
+    
+    A = A1/(4*pi);
+
+end
\ No newline at end of file

diff --git a/src/gauss.m b/src/gauss.m
new file mode 100644 (file)
index 0000000..48a1099
--- /dev/null
+++ b/src/gauss.m
@@ -0,0 +1,40 @@
+function [nodes,weights] = gauss(n,varargin)
+
+% Compute Gauss-Legendre quadrature on compact interval [a,b]
+%
+% Usage: [NODES,WEIGHTS] = gauss(N [,A,B])
+%
+% Here, N is the length of the Gaussian quadrature rule, and the
+% optional scalars A and B determine the compact interval. By default,
+% there holds A = -1 and b = +1.
+%
+% The function returns an (N x 1) column vector WEIGHTS containing the
+% quadrature weights and an (1 x N) row vector NODES containing the
+% corresponding quadrature nodes.
+%
+% Example: Suppose we aim to approximate the integral 
+%
+%    int_a^b f dx
+%
+% Assume that F is a Matlab function, which takes a column vector of
+% evaluation points X and returns a column vector Y of function values,
+% where Y(j) = F(X(j)). Then, the numerical integration reads
+%
+%    [NODES,WEIGHTS] = gauss(n,A,B);
+%    INT = WEIGHTS*F(NODES);
+%
+% Author: Dirk Praetorius - 11.03.2010
+
+
+beta = (1:n-1)./sqrt((2*(1:n-1)).^2-1);
+A = diag(beta,-1)+diag(beta,1);
+[eigenvector,nodes] = eig(A);
+[nodes,idx] = sort(diag(nodes));
+weights = 2*eigenvector(1,idx).^2;
+
+if nargin >= 3
+    a = varargin{1};
+    b = varargin{2};
+    weights = 0.5*abs(b-a)*weights;
+    nodes = 0.5*( a+b + nodes*(b-a) );
+end

diff --git a/src/surfDoubleQuad.m b/src/surfDoubleQuad.m
new file mode 100644 (file)
index 0000000..85970e2
--- /dev/null
+++ b/src/surfDoubleQuad.m
@@ -0,0 +1,31 @@
+function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
+%
+% in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
+%
+% Berechnet das Vierfachintegral über die Funktion f(x1,x2,y1,y2) mit den 
+% Grenzen [a b]x[c d]x[r t]x[u v]. Dazu wird eine Gaussquadratur über w 
+% Punkte verwendet.
+%
+% P. Schaefer
+
+    [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)
+            parfor 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
\ No newline at end of file