From: Peter Schaefer Date: Thu, 29 Mar 2012 09:22:44 +0000 (+0200) Subject: [src] Quadratur wieder aktiviert X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=2debe566af68f4abcbb45e8210e4294eeb508955;p=bacc.git [src] Quadratur wieder aktiviert [src] Quadratur wieder angepasst --- diff --git a/src/A_plots.m b/src/A_plots.m index 0b603b7..a612e6c 100644 --- 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 71657c1..038d47f 100644 --- 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 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 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 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