From 6a1941f09b8b387cbe36be8911b3663ddecad3db Mon Sep 17 00:00:00 2001 From: treecity Date: Mon, 27 Jun 2011 13:07:18 +0000 Subject: [PATCH] =?utf8?q?[src]=20SLPrecangle=20umbenannt=20[src]=20quadIn?= =?utf8?q?t=20in=20applyInt=20umbenannt=20[doc]=20beschreibungen=20in=20ei?= =?utf8?q?nigen=20Matlab-Funktionen=20hinzugef=C3=BCgt?= MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@34 26120e32-c555-405d-b3e1-1f783fb42516 --- src/mex_Fort.cpp | 2 +- src/mex_Fpar.cpp | 2 +- src/mex_G00.cpp | 2 +- src/{build_A.cpp => mex_build_A.cpp} | 14 ++-- src/mex_g0.cpp | 2 +- src/mexit.m | 10 +-- src/quadNorm.m | 1 + src/refineQuad.m | 14 ++++ src/{SLPrecangle.cpp => slpRectangle.cpp} | 6 +- src/{SLPrecangle.hpp => slpRectangle.hpp} | 4 +- src/surfDoubleQuad.m | 96 +++-------------------- src/surfQuad.m | 29 ++----- src/test_solve.m | 2 +- src/test_solveError.m | 4 +- src/test_solveS.m | 2 +- 15 files changed, 59 insertions(+), 131 deletions(-) rename src/{build_A.cpp => mex_build_A.cpp} (91%) rename src/{SLPrecangle.cpp => slpRectangle.cpp} (96%) rename src/{SLPrecangle.hpp => slpRectangle.hpp} (97%) diff --git a/src/mex_Fort.cpp b/src/mex_Fort.cpp index 7d4f894..0f4b141 100644 --- a/src/mex_Fort.cpp +++ b/src/mex_Fort.cpp @@ -4,7 +4,7 @@ #include "mex.h" #include -#include "SLPrecangle.hpp" +#include "slpRectangle.hpp" #define my_PI 3.141592653589793 #define EPS 0.02 diff --git a/src/mex_Fpar.cpp b/src/mex_Fpar.cpp index ab8ae94..c02cfb6 100644 --- a/src/mex_Fpar.cpp +++ b/src/mex_Fpar.cpp @@ -4,7 +4,7 @@ #include "mex.h" #include -#include "SLPrecangle.hpp" +#include "slpRectangle.hpp" #define my_PI 3.141592653589793 #define EPS 0.02 diff --git a/src/mex_G00.cpp b/src/mex_G00.cpp index 0bc77f3..9f9fd98 100644 --- a/src/mex_G00.cpp +++ b/src/mex_G00.cpp @@ -4,7 +4,7 @@ #include "mex.h" #include -#include "SLPrecangle.hpp" +#include "slpRectangle.hpp" #define my_PI 3.141592653589793 #define EPS 0.02 diff --git a/src/build_A.cpp b/src/mex_build_A.cpp similarity index 91% rename from src/build_A.cpp rename to src/mex_build_A.cpp index 7268af9..48a30b0 100644 --- a/src/build_A.cpp +++ b/src/mex_build_A.cpp @@ -7,7 +7,7 @@ //#include "tbb/parallel_for.h" -#include "SLPrecangle.hpp" +#include "slpRectangle.hpp" #define EPS 0.02 @@ -225,7 +225,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { // y1[0], y0[1], y3[1], d[0], d[1], d[2]); //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]); tmp - = quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]), + = apply0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]), fabs(yb[rxb]), d[rxa], d[rxb], d[rx]); // printf("%d%d|%.2f\n",j,k,tmp); @@ -234,7 +234,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { // y1[0], y0[1], y3[1], d[0], d[1], d[2]); //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]); tmp - = quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]), + = apply0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[rxa]), fabs(ya[rxb]), d[rxa], d[rxb], d[rx]); // printf("%d%d|%.2f\n",j,k,tmp); @@ -247,22 +247,22 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (rxa == rya) { tmp - = quad0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]), + = apply0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]), fabs(ya[rya]), fabs(yb[ryb]), d[rxb], d[rxa], d[rx]); } else if (rxa == ryb) { tmp - = quad0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]), + = apply0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]), fabs(yb[ryb]), fabs(ya[rya]), d[rxb], d[rxa], d[rx]); } else if (rxb == rya) { tmp - = quad0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]), + = apply0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rya]), fabs(yb[ryb]), d[rxa], d[rxb], d[rx]); } else { tmp - = quad0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]), + = apply0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[ryb]), fabs(ya[rya]), d[rxa], d[rxb], d[rx]); } diff --git a/src/mex_g0.cpp b/src/mex_g0.cpp index a93de87..9924813 100644 --- a/src/mex_g0.cpp +++ b/src/mex_g0.cpp @@ -4,7 +4,7 @@ #include "mex.h" #include -#include "SLPrecangle.hpp" +#include "slpRectangle.hpp" #define my_PI 3.141592653589793 #define EPS 0.02 diff --git a/src/mexit.m b/src/mexit.m index f552eff..ed3eaad 100644 --- a/src/mexit.m +++ b/src/mexit.m @@ -1,8 +1,8 @@ -mex mex_g0.cpp SLPrecangle.cpp -mex mex_G00.cpp SLPrecangle.cpp -mex mex_Fpar.cpp SLPrecangle.cpp -mex mex_Fort.cpp SLPrecangle.cpp +mex mex_g0.cpp slpRectangle.cpp +mex mex_G00.cpp slpRectangle.cpp +mex mex_Fpar.cpp slpRectangle.cpp +mex mex_Fort.cpp slpRectangle.cpp -mex build_A.cpp SLPrecangle.cpp +mex mex_build_A.cpp slpRectangle.cpp diff --git a/src/quadNorm.m b/src/quadNorm.m index c663ae4..18b0f85 100644 --- a/src/quadNorm.m +++ b/src/quadNorm.m @@ -6,6 +6,7 @@ function n = quadNorm(coordinates, elements,varargin) % Diese Funktion Berechnet die Orthogonalen mit Laenge 1 über alle Flächen % FLAG: % w -> Laenge entspricht Flaecheninhalt +% % P.Schaefer %% Parameterueberpruefung diff --git a/src/refineQuad.m b/src/refineQuad.m index 8502634..aa4c2ae 100644 --- a/src/refineQuad.m +++ b/src/refineQuad.m @@ -1,4 +1,18 @@ function [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type) +% +% [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type) +% +% Verfeinert die markierten Elemente mit dem entsprechenden TYP und gibt +% auch die F2S Beziehungen zurück. type muss von der Länge der Anzahl der +% Elemente entsprechen und die Einträge können 1,2,3,4 sein. +% +% der Typ zu jedem Element entspricht dabei: +% 1- keine Verfeinerung +% 2- 4 neue Elemente +% 3- 2 neue Elemente, übereinander +% 4- 2 neue Elemente, nebeneinander +% +% P. Schaefer if([1 1] == size(type)) type = repmat(type, size(elements,1),1); diff --git a/src/SLPrecangle.cpp b/src/slpRectangle.cpp similarity index 96% rename from src/SLPrecangle.cpp rename to src/slpRectangle.cpp index faefd8d..7e9a7b0 100644 --- a/src/SLPrecangle.cpp +++ b/src/slpRectangle.cpp @@ -3,7 +3,7 @@ #include #include -#include "SLPrecangle.hpp" +#include "slpRectangle.hpp" #define EPS 0.00001 @@ -134,7 +134,7 @@ double F_ort(double x1, double x2, double y2, double y3, double d1, double d2, return sol / 2.; } -double quadInt( +double applyInt( double(*f)(double, double, double, double, double, double, double), double a, double b, double c, double d, double r, double t, double u, double v, double d1, double d2, double d3) { @@ -159,7 +159,7 @@ double quadInt( */ } -double quad0Int( +double apply0Int( double(*f)(double, double, double, double, double, double, double), double b, double d, double t, double v, double d1, double d2, double d3) { diff --git a/src/SLPrecangle.hpp b/src/slpRectangle.hpp similarity index 97% rename from src/SLPrecangle.hpp rename to src/slpRectangle.hpp index d801cfe..26ebd71 100644 --- a/src/SLPrecangle.hpp +++ b/src/slpRectangle.hpp @@ -20,13 +20,13 @@ double FLO_plane(double, double, double, double, double, double, double); double FLO_perp(double, double, double, double, double, double, double); // sol = quadInt((F_par/F_ort),a,b,c,d,r,t,u,v,d1,d2,d3); -double quadInt( +double applyInt( double(*)(double, double, double, double, double, double, double), double, double, double, double, double, double, double, double, double, double, double); // sol = quad0Int((F_par/F_ort),b,d,t,v,d1,d2,d3); -double quad0Int( +double apply0Int( double(*)(double, double, double, double, double, double, double), double, double, double, double, double, double, double); diff --git a/src/surfDoubleQuad.m b/src/surfDoubleQuad.m index 86ed43c..4759b7a 100644 --- a/src/surfDoubleQuad.m +++ b/src/surfDoubleQuad.m @@ -1,51 +1,12 @@ -% function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w) -% %Vierfachintegral ueber f(x1,x2,y1,y2) auf [a b]x[c d]x[r t]x[u v] mit -% %der Simpsonregel und w Auswertungspunkten -% s = (b-a)/w; -% in = 0; -% for i = a:s:b-s -% in = in + fx2(f,i,c,d,r,t,u,v,w) + ... -% 4*fx2(f,(2*i+s)/2,c,d,r,t,u,v,w) + ... -% fx2(f,i+s,c,d,r,t,u,v,w); -% end -% in = in *s/6; -% end -% -% function in = fx2(f,x1,c,d,r,t,u,v,w) -% s = (d-c)/w; -% in = 0; -% for i = c:s:d-s -% in = in + fy1(f,x1,i,r,t,u,v,w) + ... -% 4*fy1(f,x1,(2*i+s)/2,r,t,u,v,w) + ... -% fy1(f,x1,i+s,r,t,u,v,w); -% end -% in = in *s/6; -% end -% -% function in = fy1(f,x1,x2,r,t,u,v,w) -% s = (t-r)/w; -% in = 0; -% for i = r:s:t-s -% in = in + fy2(f,x1,x2,i,u,v,w) + ... -% 4*fy2(f,x1,x2,(2*i+s)/2,u,v,w) + ... -% fy2(f,x1,x2,i+s,u,v,w); -% end -% in = in *s/6; -% end -% -% function in = fy2(f,x1,x2,y1,u,v,w) -% s = (v-u)/w; -% in = 0; -% for i = u:s:v-s -% in = in + f(x1,x2,y1,i) + ... -% 4*f(x1,x2,y1,(2*i+s)/2) + ... -% f(x1,x2,y1,i+s); -% end -% in = in *s/6; -% end - - 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); @@ -58,47 +19,12 @@ function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w) 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)); + 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); -% -% -% d1=r-a; -% d2=u-c; -% % sum(x1w)a -% x1n=x1n-d1; x2n=x2n-d2; -% % c -% % d -% -% p = -0.5; -% -% in = 0; -% -% for k=1:length(x2n) -% for l = 1:length(x1n) -% -% [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 -% % in -% -% end \ No newline at end of file +end \ No newline at end of file diff --git a/src/surfQuad.m b/src/surfQuad.m index 3111d7d..36b2289 100644 --- a/src/surfQuad.m +++ b/src/surfQuad.m @@ -1,25 +1,12 @@ -% function in = surfQuad(f,a,b,c,d,v) -% %Doppelintegral ueber f(x,y) auf [a b]x[c d] mit der Simpsonregel und v -% %Auswertungspunkten -% s = (b-a)/v; -% in = 0; -% for i = a:s:b-s -% in = in + fy(f,i,c,d,v) + 4*fy(f,(2*i+s)/2,c,d,v) + fy(f,i+s,c,d,v); -% end -% in = in *s/6; -% end -% -% function in = fy(f,x,c,d,v) -% s = (d-c)/v; -% in = 0; -% for i = c:s:d-s -% in = in + f(x,i) + 4*f(x,(2*i+s)/2) + f(x,i+s); -% end -% in = in *s/6; -% end - - function in = surfQuad(f,a,b,c,d,v) +% +% in = surfQuad(f,a,b,c,d,v) +% +% Berechnet das Doppelintegral über die Funktion f(x,y) mit den Grenzen +% [a b]x[c d]. Dazu wird eine Gaussquadratur über v Punkte verwendet. +% +% P. Schaefer + [xn xw] = gauss(v,a,b); [yn yw] = gauss(v,c,d); diff --git a/src/test_solve.m b/src/test_solve.m index 184dc97..c6245c7 100644 --- a/src/test_solve.m +++ b/src/test_solve.m @@ -2,7 +2,7 @@ N = size(elements,1); tic -A = build_A(coordinates,elements); +A = mex_build_A(coordinates,elements); t = toc; b = sqrt(sum(quadNorm(coordinates,elements,'w').^2,2)); diff --git a/src/test_solveError.m b/src/test_solveError.m index fa16e15..b9a94a5 100644 --- a/src/test_solveError.m +++ b/src/test_solveError.m @@ -16,7 +16,7 @@ while size(elements,1)