From 8ddaa7e97f08d795934a23b798f3b60060c83679 Mon Sep 17 00:00:00 2001 From: treecity Date: Mon, 31 Oct 2011 08:59:18 +0000 Subject: [PATCH] =?utf8?q?[src]=20neue=20Test's=20mit=20einem=20Element=20?= =?utf8?q?[src]=20mexTest=20->=20l=C3=A4uft=20MEX=20=C3=BCberhaupt?= 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@63 26120e32-c555-405d-b3e1-1f783fb42516 --- src/mex_build_AU.cpp | 18 ++++++++++------ src/mex_hello.cpp | 9 ++++++++ src/slpRectangle.cpp | 51 ++++++++++++++++++++++++++++---------------- src/test_calcInt1.m | 1 + src/test_calcInt2.m | 44 ++++++++++++++++++++++++++++++++++++++ src/test_calcInt3.m | 47 ++++++++++++++++++++++++++++++++++++++++ src/test_calcInt4.m | 46 +++++++++++++++++++++++++++++++++++++++ 7 files changed, 191 insertions(+), 25 deletions(-) create mode 100644 src/mex_hello.cpp create mode 100644 src/test_calcInt1.m create mode 100644 src/test_calcInt2.m create mode 100644 src/test_calcInt3.m create mode 100644 src/test_calcInt4.m diff --git a/src/mex_build_AU.cpp b/src/mex_build_AU.cpp index a57c197..3dc3aae 100644 --- a/src/mex_build_AU.cpp +++ b/src/mex_build_AU.cpp @@ -21,14 +21,14 @@ using namespace std; int dimOfVec(double* vec) { - if (fabs(vec[2]) > M_EPS) + if (vec[2] != 0) return 2; - else if (fabs(vec[1]) > M_EPS) + else if (vec[1] !=0) return 1; - else if (fabs(vec[0]) > M_EPS) + else if (vec[0] !=0) return 0; else { - cerr << "dimOfVec : (" << vec[0] << vec[1] << vec[2] + cerr << "dimOfVec : (" << vec[0] << " " << vec[1]<< " " << vec[2] << ") alle Werte 0 \n"; return -1; } @@ -105,22 +105,26 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { case 0: ctypeP = calcParIntO0; ctypeO = calcOrtIntO0; + break; case 1: ctypeP = calcParIntO1; ctypeO = calcOrtIntO1; + break; case 2: ctypeP = calcParIntO2; ctypeO = calcOrtIntO2; + break; case 3: ctypeP = calcParIntO3; ctypeO = calcOrtIntO3; + break; } //LageInformationen int rx, rxa, rxb, ry, rya, ryb; - cout << " Progress: #"; +// cout << " Progress: #"; count = 0; //Ausrechnen for (j = 0; j < em; ++j) { @@ -263,11 +267,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { A[(k * em) + j] = 1. / (4 * M_PI) * tmp; if(k!=j) A[(j * em) + k] = 1. / (4 * M_PI) * tmp; - if(count++ > ((em*(em+1))/2)/10){ +/* if(count++ > ((em*(em+1))/2)/10){ count = 0; cout << "#"; cout.flush(); - } + }*/ } diff --git a/src/mex_hello.cpp b/src/mex_hello.cpp new file mode 100644 index 0000000..7911037 --- /dev/null +++ b/src/mex_hello.cpp @@ -0,0 +1,9 @@ +#include +#include "mex.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { + std::cout << "Hello" << std::endl; + + + +} diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp index 6a29495..d984cf8 100644 --- a/src/slpRectangle.cpp +++ b/src/slpRectangle.cpp @@ -283,7 +283,6 @@ double calcParIntQY1X1(double b, double d, double t, double v, double d1, } return sol; - } double calcParIntQY1(double b, double d, double t, double v, double d1, @@ -299,9 +298,25 @@ double calcOrtIntA(double b, double d, double t, double v, double d1, } -double calcOrtIntQY1Y2(double b, double d, double t, double v, double d1, +double calcOrtIntQX1X2(double b, double d, double t, double v, double d1, double d2, double d3) { - return 0; + //Fall 2 + double sol = 0; + + for (int i = 0; i < 5; ++i) { + for (int j = 0; j < 5; ++j) { + sol += G_AY1Y2(-0.5, t + d2, v + d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1) + * b * d * gauss_w5[i] * gauss_w5[j]; + sol -= G_AY1Y2(-0.5, d2, v + d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1) + * b * d * gauss_w5[i] * gauss_w5[j]; + sol -= G_AY1Y2(-0.5, t + d2, d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1) + * b * d * gauss_w5[i] * gauss_w5[j]; + sol += G_AY1Y2(-0.5, d2, d3, d * gauss_n5[j],0, b * gauss_n5[i] - d1) + * b * d * gauss_w5[i] * gauss_w5[j]; + } + } + + return sol; } @@ -344,11 +359,11 @@ double calcOrtIntO1(double b, double d, double t, double v, double d1, double tmp = 0; tmp = b; - b = t; - t = tmp; - tmp = d; - d = v; + b = v; v = tmp; + tmp = d; + d = t; + t = tmp; d1 = -d1; d2 = -d2; @@ -389,15 +404,15 @@ double calcOrtIntO2(double b, double d, double t, double v, double d1, double d2, double d3, double eta) { //Elmente anordnen - if (b * b + d * d > t * t + v * v) { + if (d > t) { double tmp = 0; tmp = b; - b = t; - t = tmp; - tmp = d; - d = v; + b = v; v = tmp; + tmp = d; + d = t; + t = tmp; d1 = -d1; d2 = -d2; @@ -405,7 +420,7 @@ double calcOrtIntO2(double b, double d, double t, double v, double d1, } if ((b * b + d * d) * eta <= d1 * d1 + d2 * d2 + d3 * d3) { - return 0; + return calcOrtIntQX1X2(b, d, t, v, d1, d2, d3); } else { return calcOrtIntA(b, d, t, v, d1, d2, d3); } @@ -442,15 +457,15 @@ double calcOrtIntO3(double b, double d, double t, double v, double d1, double d2, double d3, double eta) { //Elmente anordnen - if (b * b + d * d > t * t + v * v) { + if (d > t) { double tmp = 0; tmp = b; - b = t; - t = tmp; - tmp = d; - d = v; + b = v; v = tmp; + tmp = d; + d = t; + t = tmp; d1 = -d1; d2 = -d2; diff --git a/src/test_calcInt1.m b/src/test_calcInt1.m new file mode 100644 index 0000000..5a9e679 --- /dev/null +++ b/src/test_calcInt1.m @@ -0,0 +1 @@ + elements=[1 2 3 4;5 6 7 8]; neigh = zeros(2,8); dat = []; h = 1; coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; 2 0 0; 2+h 0 0; 2+h h 0 ; 2 h 0]; figure(1) current = coordinates(elements(1,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'g'); hold on current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'b'); h = h/2; coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; 2 0 0; 2+h 0 0; 2+h h 0 ; 2 h 0]; current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'y'); hold off legend('Element1(h/2)','Element2(h/2)','Element2(h/4)') h = 2; for I = 1:50 h = h/2; coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; 2 0 0; 2+h 0 0; 2+h h 0 ; 2 h 0]; A0 = mex_build_AU(coordinates,elements,0,0); A2 = mex_build_AU(coordinates,elements,1,2); A1 = mex_build_AU(coordinates,elements,1,1); I dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)]; dat(dat<0)=0; end figure(2) loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4))); legend('Analytisch','Quad Element','Element vertauschen'); xlabel('Elementgroesse'); ylabel('Integral'); \ No newline at end of file diff --git a/src/test_calcInt2.m b/src/test_calcInt2.m new file mode 100644 index 0000000..9320aec --- /dev/null +++ b/src/test_calcInt2.m @@ -0,0 +1,44 @@ + + +elements=[1 2 3 4;5 6 7 8]; +neigh = zeros(2,8); + +dat = []; + +h = 1/2; +coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; 4 0 0; 4+h 0 0; 4+h 1 0 ; 4 1 0]; +figure(1) +current = coordinates(elements(1,[1:4,1])',:); +fill3(current(:,1),current(:,2),current(:,3),'g'); +hold on +current = coordinates(elements(2,[1:4,1])',:); +fill3(current(:,1),current(:,2),current(:,3),'b'); + +h = h/2; +coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; 4 0 0; 4+h 0 0; 4+h 1 0 ; 4 1 0]; +current = coordinates(elements(2,[1:4,1])',:); +fill3(current(:,1),current(:,2),current(:,3),'y'); +hold off +legend('Element1(h/2)','Element2(h/2)','Element2(h/4)') + +h = 2; +for I = 1:50 + + h = h/2; +coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; 4 0 0; 4+h 0 0; 4+h 1 0 ; 4 1 0]; + + A0 = mex_build_AU(coordinates,elements,0,0); + A2 = mex_build_AU(coordinates,elements,1,2); + A3 = mex_build_AU(coordinates,elements,1,3); + I + dat(I,1:4) = [h A0(1,2) A2(1,2) A3(1,2)]; + dat(dat<0)=0; +end + +figure(2) +loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4))); + +legend('Analytisch','Quad Element','Quad Achse'); + +xlabel('Elementgroesse'); +ylabel('Integral'); \ No newline at end of file diff --git a/src/test_calcInt3.m b/src/test_calcInt3.m new file mode 100644 index 0000000..f83514e --- /dev/null +++ b/src/test_calcInt3.m @@ -0,0 +1,47 @@ + + +elements=[1 2 3 4;5 6 7 8]; +neigh = zeros(2,8); + +dat = []; + +h = 1/2; +coordinates=[0 0 0;h 0 0; h 1 0;0 1 0; 2 0 0; 2+h 0 0; 2+h 1 0 ; 2 1 0]; +figure(1) +current = coordinates(elements(1,[1:4,1])',:); +fill3(current(:,1),current(:,2),current(:,3),'g'); +hold on +current = coordinates(elements(2,[1:4,1])',:); +fill3(current(:,1),current(:,2),current(:,3),'b'); + +h = h/2; +coordinates=[0 0 0;h 0 0; h 1 0;0 1 0; 2 0 0; 2+h 0 0; 2+h 1 0 ; 2 1 0]; +current = coordinates(elements(1,[1:4,1])',:); +fill3(current(:,1),current(:,2),current(:,3),'y'); +hold on +current = coordinates(elements(2,[1:4,1])',:); +fill3(current(:,1),current(:,2),current(:,3),'r'); +hold off +legend('Element1(h/2)','Element2(h/2)','Element1(h/4)','Element2(h/4)') + +h = 2; +for I = 1:50 + + h = h/2; +coordinates=[0 0 0;h 0 0; h 1 0;0 1 0; 2 0 0; 2+h 0 0; 2+h 1 0 ; 2 1 0]; + + A0 = mex_build_AU(coordinates,elements,0,0); + A2 = mex_build_AU(coordinates,elements,1,2); + A3 = mex_build_AU(coordinates,elements,1,3); + I + dat(I,1:4) = [h A0(1,2) A2(1,2) A3(1,2)]; + dat(dat<0)=0; +end + +figure(2) +loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4))); + +legend('Analytisch','Quad Element','Quad Achse'); + +xlabel('Elementgroesse'); +ylabel('Integral'); \ No newline at end of file diff --git a/src/test_calcInt4.m b/src/test_calcInt4.m new file mode 100644 index 0000000..c2b60d7 --- /dev/null +++ b/src/test_calcInt4.m @@ -0,0 +1,46 @@ + + +elements=[1 2 3 4;5 6 7 8]; +neigh = zeros(2,8); + +dat = []; + +h = 1/2; +coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; 2 0 0; 2 0 h; 2 h h ; 2 h 0]; +figure(1) +current = coordinates(elements(1,[1:4,1])',:); +fill3(current(:,1),current(:,2),current(:,3),'g'); +hold on +current = coordinates(elements(2,[1:4,1])',:); +fill3(current(:,1),current(:,2),current(:,3),'b'); + +h = h/2; +coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; 2 0 0; 2 0 h; 2 h h ; 2 h 0]; +current = coordinates(elements(2,[1:4,1])',:); +fill3(current(:,1),current(:,2),current(:,3),'r'); +hold off +legend('Element1(h/2)','Element2(h/2)','Element1(h/4)','Element2(h/4)') + +h = 2; +for I = 1:50 + + h = h/2; +coordinates=[0 0 0;1 0 0; 1 1 0;0 1 0; 2 0 0; 2 0 h; 2 h h ; 2 h 0]; + + A0 = mex_build_AU(coordinates,elements,0,0); + A1 = mex_build_AU(coordinates,elements,1,1); + A2 = mex_build_AU(coordinates,elements,1,2); + I + dat(I,1:4) = [h A0(1,2) A1(1,2) A2(1,2)]; + + dat(dat<0)=0; + +end + +figure(2) +loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),dat(:,1),abs(dat(:,4))); + +legend('Analytisch','Element vertauschen','Quad Element','location','southeast'); + +xlabel('Elementgroesse'); +ylabel('Integral'); \ No newline at end of file -- 2.47.3