From: treecity Date: Sat, 14 Jan 2012 21:12:08 +0000 (+0000) Subject: [src] mex_build_AU new Funktions X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=e1ca4a6df879a5234c2183f60cee0292e878ee3f;p=bacc.git [src] mex_build_AU new Funktions [src] test_... doing what they should [src] A_plot added... git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@72 26120e32-c555-405d-b3e1-1f783fb42516 --- diff --git a/.cproject b/.cproject index 24c6a1f..b3be917 100644 --- a/.cproject +++ b/.cproject @@ -6,7 +6,18 @@ - + + + + + + + + + + + + @@ -28,9 +39,7 @@ @@ -38,9 +47,7 @@ @@ -48,7 +55,9 @@ + @@ -56,11 +65,7 @@ - + @@ -97,6 +102,7 @@ + @@ -107,6 +113,7 @@ + @@ -117,6 +124,9 @@ + @@ -127,6 +137,7 @@ + diff --git a/src/A_AnIso.m b/src/A_AnIso.m index d60bbcc..74345cd 100644 --- a/src/A_AnIso.m +++ b/src/A_AnIso.m @@ -7,7 +7,7 @@ A_loadMeshF(file) global G_Ta; global G_Da; -if(length(mu)~=length(type)&&length(mu)~=1) +if(length(mu)~=max(type)-1&&length(mu)~=1) disp 'Error: Pleas set right type and mu parameters' return end @@ -16,9 +16,10 @@ for i = 1:times if(size(G_Ta,1)>2) nextF = G_Ta(end,1)*4; nextS = neville2(1:size(G_Ta,1),G_Ta(:,1)',size(G_Ta,1)+1); - calc = size(G_Ta,1)+(-1:0); + cald = size(G_Ta,1)+(-1:0); + calc = size(G_Ta,1)+(-2:0); nextTime = [nextS neville2(G_Ta(calc,1)',G_Ta(calc,2)',nextF) ... - neville2(G_Ta(calc,1)',G_Ta(calc,3)',nextF) ... + neville2(G_Ta(cald,1)',G_Ta(cald,3)',nextF) ... neville2(G_Ta(calc,1)',G_Ta(calc,4)',nextS)] %overFlowTime = neville2(G_Ta(size(G_Ta,1)+(-1:0),1)',sum(G_Ta(size(G_Ta,1)+(-1:0),2:4),2)',4500) diff --git a/src/A_Iso.m b/src/A_Iso.m index bf51d6c..8f78739 100644 --- a/src/A_Iso.m +++ b/src/A_Iso.m @@ -1,13 +1,14 @@ function A_Iso(file,times,mu,type) %A_Iso(file,times,mu,type) +format longG A_loadMeshF(file) global G_Ti; global G_Di; -if(length(mu)~=length(type)&&length(mu)~=1) +if(length(mu)~=max(type)-1&&length(mu)~=1) disp 'Error: Pleas set right type and mu parameters' return end @@ -15,8 +16,9 @@ end for i = 1:times if(size(G_Ti,1)>2) nextF = G_Ti(end,1)*4; - calc = size(G_Ti,1)+(-1:0); - nextTime = [nextF neville2(G_Ti(calc,1)',G_Ti(calc,2)',nextF) ... + cald = size(G_Ti,1)+(-1:0); + calc = size(G_Ti,1)+(-2:0); + nextTime = [nextF neville2(G_Ti(cald,1)',G_Ti(cald,2)',nextF) ... neville2(G_Ti(calc,1)',G_Ti(calc,3)',nextF)] %overFlowTime = neville2(G_Ti(size(G_Ti,1)+(-1:0),1)',sum(G_Ti(size(G_Ti,1)+(-1:0),2:4),2)',4500) diff --git a/src/A_plot.m b/src/A_plot.m new file mode 100644 index 0000000..9789e21 --- /dev/null +++ b/src/A_plot.m @@ -0,0 +1,62 @@ +function A_plot() + +type2str = ['Analytisch ' ; 'Quad Element' ; 'Quad Achse '; 'Quad Seite ']; + + +%% G_Di +global G_Di; + +[m n] = size(G_Di); + +step = round(n/3); + + +if step<1 + disp ('Error: No Data to show.') +else + +figure(2) +loglog(G_Di(:,1),G_Di(:,[3+(0:step-1)*3])) + +xlabel('Elemente'); +ylabel('Schaetzer'); +legend(type2str (G_Di(1,[2+(0:step-1)*3])',:)); + + +figure(3) +loglog(G_Di(:,1),G_Di(:,[4+(0:step-1)*3])) + +xlabel('Elemente'); +ylabel('eNorm'); +legend(type2str(G_Di(1,[2+(0:step-1)*3])',:)); +end + +%%G_Da +global G_Da; + +[m n] = size(G_Da); + +step = round(n/3); + + +if step<1 + disp ('Error: No Data to show.') +else + +figure(4) +loglog(G_Da(:,1),G_Da(:,[3+(0:step-1)*3])) + +xlabel('Elemente'); +ylabel('Schaetzer'); +legend(type2str (G_Da(1,[2+(0:step-1)*3])',:)); + + +figure(5) +loglog(G_Da(:,1),G_Da(:,[4+(0:step-1)*3])) + +xlabel('Elemente');3 +ylabel('eNorm'); +legend(type2str(G_Da(1,[2+(0:step-1)*3])',:)); +end + +end \ No newline at end of file diff --git a/src/A_stepAnIso.m b/src/A_stepAnIso.m index 4b3fa1f..35d58aa 100644 --- a/src/A_stepAnIso.m +++ b/src/A_stepAnIso.m @@ -12,12 +12,12 @@ if (isempty(G_E) || isempty(G_C) || isempty(G_N)) return end -if(length(mu)~=length(type)&&length(mu)~=1) +if(length(mu)~=max(type)-1&&length(mu)~=1) disp 'Error: Pleas set right type and mu parameters' return end if(length(mu)==1) - mu = repmat(mu,length(type),1); + mu = repmat(mu,max(type)-1,1); end time = zeros(1,3); @@ -31,7 +31,7 @@ time = zeros(1,3); tic dataAniso = size(G_E,1); for i = 1:length(type) - A_fine = mex_build_AU(coordinates_fine,elements_fine,mu(i),type(i)); + A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type(i)); x_fine = A_fine\b; diff --git a/src/A_stepIso.m b/src/A_stepIso.m index 41ddd61..1abcfad 100644 --- a/src/A_stepIso.m +++ b/src/A_stepIso.m @@ -12,12 +12,12 @@ if (isempty(G_E) || isempty(G_C)) return end -if(length(mu)~=length(type)&&length(mu)~=1) +if(length(mu)~=max(type)-1&&length(mu)~=1) disp 'Error: Pleas set right type and mu parameters' return end if(length(mu)==1) - mu = repmat(mu,length(type),1); + mu = repmat(mu,max(type)-1,1); end time = zeros(1,2); @@ -30,7 +30,7 @@ time = zeros(1,2); tic dataIso = size(G_E,1); for i = 1:length(type) - A_fine = mex_build_AU(coordinates_fine,elements_fine,mu(i),type(i)); + A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type(i)); [r c] = find(isnan(A_fine)~=isinf(A_fine)); if(~isempty(r)) figure(9) diff --git a/src/mex_Int.cpp b/src/mex_Int.cpp new file mode 100644 index 0000000..5b16086 --- /dev/null +++ b/src/mex_Int.cpp @@ -0,0 +1,297 @@ +/* + * Art der Berechnung durch Parameter bestimmt... + * Analytisch/Semianalytisch/... + * + * + */ + +#include +#include +#include +#include "mex.h" +#include + +//#include "gauss.hpp" + +#define M_EPS 10e-8 + +//#include "tbb/parallel_for.h" + + +#include "slpRectangle.hpp" + +using namespace std; +using namespace slpR; + +int dimOfVec(double* vec) { + if (vec[2] != 0) + return 2; + else if (vec[1] !=0) + return 1; + else if (vec[0] !=0) + return 0; + else { + cerr << "dimOfVec : (" << vec[0] << " " << vec[1]<< " " << vec[2] + << ") alle Werte 0 \n"; + return -1; + } + +} + +inline int dimOfThird(int a, int b) { + return ((-(a + b) % 3) + 3) % 3; +} + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { + +// initeQuad(); +// cout << Quad[1].nod[0]; + + int i, j, k; //Schleifenindizes + double tmp; //Zwischenspeicherung der Einzelnen Werte + int count; + + //Sicherheitsabfragen zu Datengroessen + if ((nrhs != 4)) + mexErrMsgTxt( + "expected (coordinates(Nx3),elements(Mx4),mu(double),type(int))"); + if (nlhs > 1) + mexErrMsgTxt("has only one output argument"); + + int cm = mxGetM(prhs[0]); + int cn = mxGetN(prhs[0]); + if (cn != 3) + mexErrMsgTxt("expected coordinates (Nx3)"); + cout << " Coordinaten:" << cm << endl; + + int em = mxGetM(prhs[1]); + int en = mxGetN(prhs[1]); + if (en != 4) + mexErrMsgTxt("expected elements (Mx4)"); + cout << " Elemente:" << em << endl; + + //Auslesen der Parameter + + plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL); + double * A = mxGetPr(plhs[0]); + double * C = mxGetPr(prhs[0]); + double * E = mxGetPr(prhs[1]); + + double mu = *(mxGetPr(prhs[2])); + int type = (int) *(mxGetPr(prhs[3])); + + //Initialisieren der Hilfsvektoren + + double * x0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + + double * y0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + + double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + + // Welche Funktion soll verwendet werden + double(*ctypeP)(double, double, double, double, double, double, double, + double); + double(*ctypeO)(double, double, double, double, double, double, double, + double); + + //Art der Berechnung bestimmen + cout << " Type:" << type << endl; + switch (type) { + 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: #"; + count = 0; + //Ausrechnen + for (j = 0; j < em; ++j) { + x0[0] = C[(int) E[j] - 1]; + x0[1] = C[cm + (int) E[j] - 1]; + x0[2] = C[2 * cm + (int) E[j] - 1]; + + x1[0] = C[(int) E[em + j] - 1]; + x1[1] = C[cm + (int) E[em + j] - 1]; + x1[2] = C[2 * cm + (int) E[em + j] - 1]; + + x2[0] = C[(int) E[2 * em + j] - 1]; + x2[1] = C[cm + (int) E[2 * em + j] - 1]; + x2[2] = C[2 * cm + (int) E[2 * em + j] - 1]; + + x3[0] = C[(int) E[3 * em + j] - 1]; + x3[1] = C[cm + (int) E[3 * em + j] - 1]; + x3[2] = C[2 * cm + (int) E[3 * em + j] - 1]; + + for (i = 0; i < 3; ++i) + xa[i] = x1[i] - x0[i]; + + for (i = 0; i < 3; ++i) + xb[i] = x3[i] - x0[i]; + + // Lageeigenschaften des Flaechenstuecks + + rxa = dimOfVec(xa); + rxb = dimOfVec(xb); + rx = dimOfThird(rxa, rxb); + + //kleinste Ecke finden und fuer \delta verwenden + + if (xa[rxa] > 0) { + if (xb[rxb] > 0) { + for (i = 0; i < 3; ++i) + dt[i] = -x0[i]; + } else { + for (i = 0; i < 3; ++i) + dt[i] = -x3[i]; + } + } else { + if (xb[rxb] > 0) { + for (i = 0; i < 3; ++i) + dt[i] = -x1[i]; + } else { + for (i = 0; i < 3; ++i) + dt[i] = -x2[i]; + } + } + + for (k = i+1; k < em; ++k) { + y0[0] = C[(int) E[k] - 1]; + y0[1] = C[cm + (int) E[k] - 1]; + y0[2] = C[2 * cm + (int) E[k] - 1]; + + y1[0] = C[(int) E[em + k] - 1]; + y1[1] = C[cm + (int) E[em + k] - 1]; + y1[2] = C[2 * cm + (int) E[em + k] - 1]; + + y2[0] = C[(int) E[2 * em + k] - 1]; + y2[1] = C[cm + (int) E[2 * em + k] - 1]; + y2[2] = C[2 * cm + (int) E[2 * em + k] - 1]; + + y3[0] = C[(int) E[3 * em + k] - 1]; + y3[1] = C[cm + (int) E[3 * em + k] - 1]; + y3[2] = C[2 * cm + (int) E[3 * em + k] - 1]; + + for (i = 0; i < 3; ++i) + ya[i] = y1[i] - y0[i]; + + for (i = 0; i < 3; ++i) + yb[i] = y3[i] - y0[i]; + + // Lageeigenschaften des Flaechenstuecks + + rya = dimOfVec(ya); + ryb = dimOfVec(yb); + ry = dimOfThird(rya, ryb); + + //kleinste Ecke finden und fuer \delta verwenden + + if (ya[rya] > 0) { + if (yb[ryb] > 0) { + for (i = 0; i < 3; ++i) + d[i] = dt[i] + y0[i]; + } else { + for (i = 0; i < 3; ++i) + d[i] = dt[i] + y3[i]; + } + } else { + if (yb[ryb] > 0) { + for (i = 0; i < 3; ++i) + d[i] = dt[i] + y1[i]; + } else { + for (i = 0; i < 3; ++i) + d[i] = dt[i] + y2[i]; + } + } + + if (rx == ry) { + if (rxa == rya) { + tmp + = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), + fabs(ya[rxa]), fabs(yb[rxb]), d[rxa], + d[rxb], d[rx], mu); + + } else { + tmp + = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), + fabs(yb[rxa]), fabs(ya[rxb]), d[rxa], + d[rxb], d[rx], mu); + } + + } else { + + if (rxa == rya) { + tmp + = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), + fabs(ya[rya]), fabs(yb[ryb]), d[rxb], + d[rxa], d[rx], mu); + } else if (rxa == ryb) { + tmp + = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), + fabs(yb[ryb]), fabs(ya[rya]), d[rxb], + d[rxa], d[rx], mu); + } else if (rxb == rya) { + tmp + = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), + fabs(ya[rya]), fabs(yb[ryb]), d[rxa], + d[rxb], d[rx], mu); + } else { + tmp + = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), + fabs(yb[ryb]), fabs(ya[rya]), d[rxa], + d[rxb], d[rx], mu); + } + + } + 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){ + count = 0; + cout << "#"; + cout.flush(); + }*/ + + } + + } + cout << endl; + //Rueckgabe (eventuell zurueck schreiben) + //mxFree(x0); + //mxFree(x1); + //mxFree(x3); + //mxFree(y0); + //mxFree(y1); + //mxFree(y3); + //mxFree(d); + //mxFree(dt); + + return; +} diff --git a/src/mex_build_AU.cpp b/src/mex_build_AU.cpp index c6b866f..7f54061 100644 --- a/src/mex_build_AU.cpp +++ b/src/mex_build_AU.cpp @@ -77,8 +77,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double * C = mxGetPr(prhs[0]); double * E = mxGetPr(prhs[1]); - double mu = *(mxGetPr(prhs[2])); + int type = (int) *(mxGetPr(prhs[3])); + double * mu = mxGetPr(prhs[2]); //Initialisieren der Hilfsvektoren @@ -101,30 +102,25 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { // Welche Funktion soll verwendet werden double(*ctypeP)(double, double, double, double, double, double, double, - double); + double*); double(*ctypeO)(double, double, double, double, double, double, double, - double); + double*); //Art der Berechnung bestimmen cout << " Type:" << type << endl; switch (type) { - case 0: - ctypeP = calcParIntO0; - ctypeO = calcOrtIntO0; - break; - case 1: - ctypeP = calcParIntO1; - ctypeO = calcOrtIntO1; + default: + ctypeP = cParO1; + ctypeO = NULL; break; case 2: - ctypeP = calcParIntO2; - ctypeO = calcOrtIntO2; + ctypeP = cParO2; + ctypeO = NULL; break; case 3: - ctypeP = calcParIntO3; - ctypeO = calcOrtIntO3; + ctypeP = cParO3; + ctypeO = NULL; break; - } //LageInformationen diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp index 2ca81a3..ccdfc86 100644 --- a/src/slpRectangle.cpp +++ b/src/slpRectangle.cpp @@ -147,7 +147,7 @@ double slpR::G_AY1Y2(double p, double y1, double y2, double x1, double x2, doubl } else { sol = NAN; cout << "warning in G_AY1Y2: no case for p=" << p - << " defined. Possible: [-1.5,-0.5,0.5,...]\n"; + << " defined. Possible: [-1.5,-0.5,0.5,...]" << endl; } return sol; @@ -157,8 +157,6 @@ double slpR::Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2, double l) { double sol = 0; - //sol += y2*G_AY1Y2(y1,x2,) - return sol; } @@ -334,6 +332,15 @@ double slpR::calcParIntQY1X1(double b, double d, double t, double v, double d1, double slpR::calcParIntQY1(double b, double d, double t, double v, double d1, double d2, double d3) { //Fall 4 + + + + + + double sol = 0; + + + return 0; ///ACHTUNG noch Falsch } @@ -438,11 +445,12 @@ double slpR::calcParIntO2(double b, double d, double t, double v, double d1, d3 = -d3; } - if ((b * b + d * d) * eta < d1 * d1 + d2 * d2 + d3 * d3) { +// if ((b * b + d * d) * eta < d1 * d1 + d2 * d2 + d3 * d3) { return calcParIntQX1X2(b, d, t, v, d1, d2, d3); - } else { - return calcParIntA(b, d, t, v, d1, d2, d3); - } +// } else { +// cout << "2error"; +// return calcParIntA(b, d, t, v, d1, d2, d3); +// } } @@ -468,6 +476,7 @@ double slpR::calcOrtIntO2(double b, double d, double t, double v, double d1, if ((b * b + d * d) * eta <= d1 * d1 + d2 * d2 + d3 * d3) { return calcOrtIntQX1X2(b, d, t, v, d1, d2, d3); } else { + cout << "2error"; return calcOrtIntA(b, d, t, v, d1, d2, d3); } } @@ -491,11 +500,12 @@ double slpR::calcParIntO3(double b, double d, double t, double v, double d1, d2 = tmp; } - if ((min(b, t) * eta <= d1) && (min(d, v) * eta <= d2) ) { +// if ((min(b, t) * eta <= d1) && (min(d, v) * eta <= d2) ) { return calcParIntQY1X1(b, d, t, v, d1, d2, d3); - } else { - return calcParIntA(b, d, t, v, d1, d2, d3); - } +// } else { +// cout << "3error"; +// return calcParIntA(b, d, t, v, d1, d2, d3); +// } } @@ -508,3 +518,95 @@ double slpR::calcOrtIntO3(double b, double d, double t, double v, double d1, return calcOrtIntA(b, d, t, v, d1, d2, d3); } } + + + +double slpR::cParO1(double b, double d, double t, double v, double d1, + double d2, double d3, double* eta) { + return calcParIntA(b, d, t, v, d1, d2, d3);; +} + + +double slpR::cParO2(double b, double d, double t, double v, double d1, + double d2, double d3, double* eta) { + + double tmp = 0; + + // kleine Seite nach vorn + if(b>t){ + tmp = b; + b = t; + t = tmp; + d1=-d1; + } + + // kleine Seite nach vorn + if(d>v){ + tmp = d; + d = v; + v = tmp; + d2=-d2; + } + + if ((b * b + d * d) * eta[0] < d1 * d1 + d2 * d2 + d3 * d3) { +// cout << "E"; + return calcParIntQX1X2(b, d, t, v, d1, d2, d3); + } else { + return calcParIntA(b, d, t, v, d1, d2, d3); + } + + return 0; +} + +double slpR::cParO3(double b, double d, double t, double v, double d1, + double d2, double d3, double* eta) { + + double tmp = 0; + + // kleine Seite nach vorn + if(b>t){ + tmp = b; + b = t; + t = tmp; + d1=-d1; + } + + // kleine Seite nach vorn + if(d>v){ + tmp = d; + d = v; + v = tmp; + d2=-d2; + } + + if ((b * b + d * d) * eta[0] < d1 * d1 + d2 * d2 + d3 * d3) { + + //kleinere Achse nach vorn + if (b * b + t * t > d * d + v * v) { + double tmp = 0; + + tmp = b; + b = d; + d = tmp; + tmp = t; + t = v; + v = tmp; + + tmp = d1; + d1 = d2; + d2 = tmp; + } + + if (min(b, t) * eta[1] <= d1){ +// cout << "A"; + return calcParIntQY1X1(b, d, t, v, d1, d2, d3); + }else{ +// cout << "E"; + return calcParIntQX1X2(b, d, t, v, d1, d2, d3); + } + } else { + return calcParIntA(b, d, t, v, d1, d2, d3); + } + + return 0; +} diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp index 801b9a6..6b56825 100644 --- a/src/slpRectangle.hpp +++ b/src/slpRectangle.hpp @@ -67,5 +67,17 @@ namespace slpR { double); + //Automatische entscheidung, welche Quad verwendet werden soll + // sol = cParIntX(b,d,t,v,d1,d2,d3,eta); + + //Voll Analytisch + double cParO1(double, double, double, double, double, double, double, + double*); + double cParO2(double, double, double, double, double, double, double, + double*); + double cParO3(double, double, double, double, double, double, double, + double*); + + } #endif diff --git a/src/test_calcInt1.m b/src/test_calcInt1.m index 2945b61..3dd4861 100644 --- a/src/test_calcInt1.m +++ b/src/test_calcInt1.m @@ -1 +1 @@ - coo = @(h,diff)[0 0 0;1 0 0; 1 1 0;0 1 0; 0 0 0 ; 0+h 0 0; 0+h h 0 ; 0 h 0]+[zeros(4,3);repmat(diff,4,1)]; elements=[1 2 3 4;5 6 7 8]; neigh = zeros(2,8); dat = []; diff = [2 0 0]; %% Laage ??bersicht figure(1) h = 1; coordinates=coo(h,diff) 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=coo(h,diff); current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'y'); hold off legend('Element1(h)','Element2(h)','Element2(h/2)'); title('Laage der Elemente'); %% Integrale bei kleiner werdenden Elementen h = 2; for I = 1:50 h = h/2; coordinates=coo(h,diff); 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','location','best'); xlabel('Elementgroesse (k??rzeste Seite)'); ylabel('Integral'); title('Integral bei kleiner werdenden Element'); %% Netzverfeinerung mit schlechtem Netz h = (1/2)^25; coordinates=coo(h,diff); A_loadMesh(coordinates,elements,neigh); datA=[]; for I = 1:4 datA(I,:) = A_stepIso(1,[0 3 2]); % datA(I,:) = A_stepAniso(1,[0 3 2],1,0); I end figure(3) loglog(datA(2:end,1),datA(2:end,3),datA(2:end,1),datA(2:end,9),datA(2:end,1),datA(2:end,6)) legend('Analytisch','Quad Element','Quad Achse','location','best'); xlabel('Anzahl der Elemente'); ylabel('mu Schaetzer'); title('mu Schaetzer mit "schlechtem" Startnetz'); figure(4) loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7)) legend('Analytisch','Quad Element','Quad Achse','location','best'); xlabel('Anzahl der Elemente'); ylabel('EnergieNorm ^2 '); title('EnergieNorm ^2 mit "schlechtem" Startnetz'); datA \ No newline at end of file + coo = @(h,diff)[0 0 0;1 0 0; 1 1 0;0 1 0; 0 0 0 ; 0+h 0 0; 0+h h 0 ; 0 h 0]+[zeros(4,3);repmat(diff,4,1)]; elements=[1 2 3 4;5 6 7 8]; neigh = zeros(2,8); dat = []; diff = [2 0 0]; %% Laage Uebersicht figure(1) h = 1; coordinates=coo(h,diff); 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=coo(h,diff); current = coordinates(elements(2,[1:4,1])',:); fill3(current(:,1),current(:,2),current(:,3),'y'); hold off legend('Element1(h)','Element2(h)','Element2(h/2)'); title('Laage der Elemente'); %% Integrale bei kleiner werdenden Elementen h = 2; for I = 1:50 h = h/2; coordinates=coo(h,diff); 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','location','best'); xlabel('Elementgroesse (k??rzeste Seite)'); ylabel('Integral'); title('Integral bei kleiner werdenden Element'); %% Netzverfeinerung mit schlechtem Netz % h = (1/2)^25; % coordinates=coo(h,diff); % % A_loadMesh(coordinates,elements,neigh); % datA=[]; % for I = 1:4 % datA(I,:) = A_stepIso(1,[0 3 2]); % % datA(I,:) = A_stepAniso(1,[0 3 2],1,0); % I % end % figure(3) % loglog(datA(2:end,1),datA(2:end,3),datA(2:end,1),datA(2:end,9),datA(2:end,1),datA(2:end,6)) % % legend('Analytisch','Quad Element','Quad Achse','location','best'); % xlabel('Anzahl der Elemente'); % ylabel('mu Schaetzer'); % title('mu Schaetzer mit "schlechtem" Startnetz'); % % figure(4) % loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7)) % % legend('Analytisch','Quad Element','Quad Achse','location','best'); % xlabel('Anzahl der Elemente'); % ylabel('EnergieNorm ^2 '); % title('EnergieNorm ^2 mit "schlechtem" Startnetz'); % % datA \ No newline at end of file diff --git a/src/test_calcInt2.m b/src/test_calcInt2.m index 6087e4b..b7300f9 100644 --- a/src/test_calcInt2.m +++ b/src/test_calcInt2.m @@ -35,9 +35,9 @@ for I = 1:50 A0 = mex_build_AU(coordinates,elements,0,0); A2 = mex_build_AU(coordinates,elements,1,2); - A1 = mex_build_AU(coordinates,elements,1,3); + A3 = mex_build_AU(coordinates,elements,1,3); I - dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)]; + dat(I,1:4) = [h A0(1,2) A2(1,2) A3(1,2)]; dat(dat<0)=0; end @@ -51,30 +51,30 @@ title('Integral bei kleiner werdenden Element'); %% Netzverfeinerung mit schlechtem Netz -h = (1/2)^25; -coordinates=coo(h,diff); - -A_loadMesh(coordinates,elements,neigh); -datA=[]; -for I = 1:4 - datA(I,:) = A_stepIso(1,[0 3 2]); -% datA(I,:) = A_stepAniso(1,[0 1 3 2],1,0); - I -end -figure(3) -loglog(datA(2:end,1),datA(2:end,3),datA(2:end,1),datA(2:end,9),datA(2:end,1),datA(2:end,6)) - -legend('Analytisch','Quad Element','Quad Achse','location','best'); -xlabel('Anzahl der Elemente'); -ylabel('mu Schaetzer'); -title('mu Schaetzer mit "schlechtem" Startnetz'); - -figure(4) -loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7)) - -legend('Analytisch','Quad Element','Quad Achse','location','best'); -xlabel('Anzahl der Elemente'); -ylabel('EnergieNorm ^2 '); -title('EnergieNorm ^2 mit "schlechtem" Startnetz'); - -datA +% h = (1/2)^25; +% coordinates=coo(h,diff); +% +% A_loadMesh(coordinates,elements,neigh); +% datA=[]; +% for I = 1:4 +% datA(I,:) = A_stepIso(1,[0 3 2]); +% % datA(I,:) = A_stepAniso(1,[0 1 3 2],1,0); +% I +% end +% figure(3) +% loglog(datA(2:end,1),datA(2:end,3),datA(2:end,1),datA(2:end,9),datA(2:end,1),datA(2:end,6)) +% +% legend('Analytisch','Quad Element','Quad Achse','location','best'); +% xlabel('Anzahl der Elemente'); +% ylabel('mu Schaetzer'); +% title('mu Schaetzer mit "schlechtem" Startnetz'); +% +% figure(4) +% loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7)) +% +% legend('Analytisch','Quad Element','Quad Achse','location','best'); +% xlabel('Anzahl der Elemente'); +% ylabel('EnergieNorm ^2 '); +% title('EnergieNorm ^2 mit "schlechtem" Startnetz'); +% +% datA diff --git a/src/test_calcInt3.m b/src/test_calcInt3.m index 359ad3d..97d7a8f 100644 --- a/src/test_calcInt3.m +++ b/src/test_calcInt3.m @@ -37,14 +37,14 @@ for I = 1:50 A0 = mex_build_AU(coordinates,elements,0,0); A2 = mex_build_AU(coordinates,elements,1,2); - A1 = mex_build_AU(coordinates,elements,1,3); + A3 = mex_build_AU(coordinates,elements,1,3); I - dat(I,1:4) = [h A0(1,2) A2(1,2) A1(1,2)]; + 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))); +loglog(dat(:,1),abs(dat(:,2)),dat(:,1),abs(dat(:,3)),'-.',dat(:,1),abs(dat(:,4)),'--'); legend('Analytisch','Quad Element','Quad Achse','location','best'); xlabel('Elementgroesse (k??rzeste Seite)'); @@ -52,31 +52,31 @@ ylabel('Integral'); title('Integral bei kleiner werdenden Element'); %% Netzverfeinerung mit schlechtem Netz - -h = (1/2)^25; -coordinates=coo(h,diff); - -A_loadMesh(coordinates,elements,neigh); -datA=[]; -for I = 1:4 - datA(I,:) = A_stepIso(1,[0 3 2]); -% datA(I,:) = A_stepAniso(1,[0 3 2],1,0); - I -end -figure(3) -loglog(datA(2:end,1),datA(2:end,3),datA(2:end,1),datA(2:end,9),datA(2:end,1),datA(2:end,6)) - -legend('Analytisch','Quad Element','Quad Achse','location','best'); -xlabel('Anzahl der Elemente'); -ylabel('mu Schaetzer'); -title('mu Schaetzer mit "schlechtem" Startnetz'); - -figure(4) -loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7)) - -legend('Analytisch','Quad Element','Quad Achse','location','best'); -xlabel('Anzahl der Elemente'); -ylabel('EnergieNorm ^2 '); -title('EnergieNorm ^2 mit "schlechtem" Startnetz'); - -datA +% +% h = (1/2)^25; +% coordinates=coo(h,diff); +% +% A_loadMesh(coordinates,elements,neigh); +% datA=[]; +% for I = 1:4 +% datA(I,:) = A_stepIso(1,[0 3 2]); +% % datA(I,:) = A_stepAniso(1,[0 3 2],1,0); +% I +% end +% figure(3) +% loglog(datA(2:end,1),datA(2:end,3),datA(2:end,1),datA(2:end,9),datA(2:end,1),datA(2:end,6)) +% +% legend('Analytisch','Quad Element','Quad Achse','location','best'); +% xlabel('Anzahl der Elemente'); +% ylabel('mu Schaetzer'); +% title('mu Schaetzer mit "schlechtem" Startnetz'); +% +% figure(4) +% loglog(datA(:,1),datA(:,4),datA(:,1),datA(:,10),datA(:,1),datA(:,7)) +% +% legend('Analytisch','Quad Element','Quad Achse','location','best'); +% xlabel('Anzahl der Elemente'); +% ylabel('EnergieNorm ^2 '); +% title('EnergieNorm ^2 mit "schlechtem" Startnetz'); +% +% datA