--- /dev/null
+function A_loadMesh(coordinates,elements)
+% Laed ein Mesh
+%
+% A_loadMesh(coordinates,elements)
+
+ global G_C;
+ global G_E;
+ G_C = coordinates;
+ G_E = elements;
+end
+
--- /dev/null
+function A_setParam(paras)
+%setzt die Parameter fuer die Durchlaeufe
+
+global G_para;
+G_para = paras;
+
+end
+
--- /dev/null
+function dataAniso = A_stepAnIso(mu,type)
+
+global G_E;
+global G_C;
+
+if (isempty(G_E) || isempty(G_C))
+ disp 'Error: Please use A_loadMesh first'
+ return
+end
+
+ [coordinates_fine,elements_fine, f2s]=refineQuad(G_C,G_E,2*ones(1,size(G_E,1)));
+
+ A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type);
+
+ b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+
+ x_fine = A_fine\b;
+
+ xe_fine = x_fine'*A_fine*x_fine;
+
+ ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
+
+ marked = mark(x_fine(f2s)',ind,0.5,0.5);
+
+ dataAniso = [size(G_E,1) sqrt(sum(ind)) xe_fine];
+ [G_C, G_E] = refineQuad(G_C,G_E,marked);
+
+
+end
+
--- /dev/null
+function dataIso = A_stepIso(mu,type)
+
+global G_E;
+global G_C;
+
+if (isempty(G_E) || isempty(G_C))
+ disp 'Error: Please use A_loadMesh first'
+ return
+end
+
+ [coordinates_fine,elements_fine, f2s]=refineQuad(G_C,G_E,2*ones(1,size(G_E,1)));
+
+ A_fine = mex_build_AU(coordinates_fine,elements_fine,mu,type);
+
+ b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+
+ x_fine = A_fine\b;
+
+ xe_fine = x_fine'*A_fine*x_fine;
+
+ ind = computeEstSlpMuTilde(x_fine,G_C,G_E,f2s);
+
+ dataIso = [size(G_E,1) sqrt(sum(ind)) xe_fine];
+ [G_C, G_E] = refineQuad(G_C,G_E,2*ones(1,size(G_E,1)));
+
+
+end
+
+++ /dev/null
-function A = build_A2(coordinates,elements)
- N = size(elements,1);
-
- A1 = zeros(N);
-
- % 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);
-
- for j = 1:N
- for k = 1:N
- ej = coordinates(elements(j,:)',:);
- ek = coordinates(elements(k,:)',:);
-
- d = ek(1,:) - ej(1,:);
-
- ej = ej - repmat(ej(1,:),3,1);
- ek = ek - repmat(ek(1,:),3,1);
-
-% d = ones(1,3);
-
-% 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(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
-
- A = A1/(4*pi);
-
-end
\ No newline at end of file
+++ /dev/null
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.hpp"
-
-#define my_PI 3.141592653589793
-#define EPS 0.02
-
-using namespace std;
-
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
- //sicherheitsabfragen zu Datengroessen
- if (nrhs != 7)
- mexErrMsgTxt("expected (x1,x2,y2,y3,d1,d2,d3)");
- if (nlhs > 2)
- mexErrMsgTxt("has maximal two output arguments");
-
-
- plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
- *mxGetPr(plhs[0]) = F_ort(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
-}
-
-
-
+++ /dev/null
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.hpp"
-
-#define my_PI 3.141592653589793
-#define EPS 0.02
-
-using namespace std;
-
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
- //sicherheitsabfragen zu Datengroessen
- if (nrhs != 7)
- mexErrMsgTxt("expected (x1,x2,y1,y2,d1,d2,d3)");
- if (nlhs > 2)
- mexErrMsgTxt("has maximal two output arguments");
-
-
- plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
- *mxGetPr(plhs[0]) = F_par(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
-}
-
-
-
+++ /dev/null
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.hpp"
-
-#define my_PI 3.141592653589793
-#define EPS 0.02
-
-using namespace std;
-
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
- //sicherheitsabfragen zu Datengroessen
- if (nrhs != 6)
- mexErrMsgTxt("expected (p, y1, y2, x1, x2, l)");
- if (nlhs > 2)
- mexErrMsgTxt("has maximal two output arguments");
-
-
- plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
- plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
- *mxGetPr(plhs[0]) = G_AY1Y2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
- *mxGetPr(plhs[1]) = G_QY1Y2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
-
-}
-
-
-
+++ /dev/null
-/*\r
- * voll Analytisch\r
- * ohne Optimierungen\r
- *\r
- *\r
- */\r
-\r
-\r
-#include <iostream>\r
-#include <cmath>\r
-#include <cassert>\r
-#include "mex.h"\r
-#include <stdlib.h>\r
-\r
-//#include "tbb/parallel_for.h"\r
-\r
-\r
-#include "slpRectangle.hpp"\r
-\r
-#define EPS 0.02\r
-\r
-using namespace std;\r
-\r
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
-\r
- int i, j, k;\r
-\r
- //sicherheitsabfragen zu Datengroessen\r
- if (nrhs != 2)\r
- mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx4))");\r
- if (nlhs > 1)\r
- mexErrMsgTxt("has only one output argument");\r
-\r
- int cm = mxGetM(prhs[0]);\r
- int cn = mxGetN(prhs[0]);\r
-\r
- if (cn != 3)\r
- mexErrMsgTxt("expected coordinates (Nx3)");\r
- int em = mxGetM(prhs[1]);\r
- int en = mxGetN(prhs[1]);\r
- if (en != 4)\r
- mexErrMsgTxt("expected elements (Mx4)");\r
- //Vorbereitung der Daten\r
-\r
- plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);\r
- double * A = mxGetPr(plhs[0]);\r
- double * C = mxGetPr(prhs[0]);\r
- double * E = mxGetPr(prhs[1]);\r
-\r
- double * x0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double * y0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * yn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double tmp;\r
-\r
- //LageInformationen\r
- int rx, rxa, rxb, ry, rya, ryb;\r
-\r
- //Ausrechnen\r
- for (j = 0; j < em; ++j) {\r
- x0[0] = C[(int) E[j] - 1];\r
- x0[1] = C[cm + (int) E[j] - 1];\r
- x0[2] = C[2 * cm + (int) E[j] - 1];\r
-\r
- x1[0] = C[(int) E[em + j] - 1];\r
- x1[1] = C[cm + (int) E[em + j] - 1];\r
- x1[2] = C[2 * cm + (int) E[em + j] - 1];\r
-\r
- x2[0] = C[(int) E[2 * em + j] - 1];\r
- x2[1] = C[cm + (int) E[2 * em + j] - 1];\r
- x2[2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
-\r
- x3[0] = C[(int) E[3 * em + j] - 1];\r
- x3[1] = C[cm + (int) E[3 * em + j] - 1];\r
- x3[2] = C[2 * cm + (int) E[3 * em + j] - 1];\r
-\r
- for (i = 0; i < 3; ++i)\r
- xa[i] = x1[i] - x0[i];\r
-\r
- for (i = 0; i < 3; ++i)\r
- xb[i] = x3[i] - x0[i];\r
-\r
- //Kreuzprodukt\r
- xn[0] = (xa[1]) * (xb[2]) - (xa[2]) * (xb[1]);\r
- xn[1] = (xa[2]) * (xb[0]) - (xa[0]) * (xb[2]);\r
- xn[2] = (xa[0]) * (xb[1]) - (xa[1]) * (xb[0]);\r
-\r
- // Lageeigenschaften des Flaechenstuecks\r
- //if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])\r
- if (xn[2] != 0)\r
- rx = 2;\r
- else if (xn[1] != 0)\r
- rx = 1;\r
- else\r
- rx = 0;\r
-\r
- if (xa[2] != 0)\r
- rxa = 2;\r
- else if (xa[1] != 0)\r
- rxa = 1;\r
- else\r
- rxa = 0;\r
-\r
- if (xb[2] != 0)\r
- rxb = 2;\r
- else if (xb[1] != 0)\r
- rxb = 1;\r
- else\r
- rxb = 0;\r
-\r
- //kleinste Ecke finden und fuer \delta verwenden\r
-\r
- if (xa[rxa] > 0) {\r
- if (xb[rxb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x0[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x3[i];\r
- }\r
- } else {\r
- if (xb[rxb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x1[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x2[i];\r
- //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
- }\r
- }\r
-\r
- //for (i=0;i<3;++i)\r
- // dt[i] = 0;\r
-\r
-\r
- // printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);\r
-\r
- for (k = 0; k < em; ++k) {\r
- y0[0] = C[(int) E[k] - 1];\r
- y0[1] = C[cm + (int) E[k] - 1];\r
- y0[2] = C[2 * cm + (int) E[k] - 1];\r
-\r
- y1[0] = C[(int) E[em + k] - 1];\r
- y1[1] = C[cm + (int) E[em + k] - 1];\r
- y1[2] = C[2 * cm + (int) E[em + k] - 1];\r
-\r
- y2[0] = C[(int) E[2 * em + k] - 1];\r
- y2[1] = C[cm + (int) E[2 * em + k] - 1];\r
- y2[2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
-\r
- y3[0] = C[(int) E[3 * em + k] - 1];\r
- y3[1] = C[cm + (int) E[3 * em + k] - 1];\r
- y3[2] = C[2 * cm + (int) E[3 * em + k] - 1];\r
-\r
- for (i = 0; i < 3; ++i)\r
- ya[i] = y1[i] - y0[i];\r
-\r
- for (i = 0; i < 3; ++i)\r
- yb[i] = y3[i] - y0[i];\r
-\r
- yn[0] = (ya[1]) * (yb[2]) - (ya[2]) * (yb[1]);\r
- yn[1] = (ya[2]) * (yb[0]) - (ya[0]) * (yb[2]);\r
- yn[2] = (ya[0]) * (yb[1]) - (ya[1]) * (yb[0]);\r
-\r
- if (yn[2] != 0)\r
- ry = 2;\r
- else if (yn[1] != 0)\r
- ry = 1;\r
- else\r
- ry = 0;\r
-\r
- if (ya[2] != 0)\r
- rya = 2;\r
- else if (ya[1] != 0)\r
- rya = 1;\r
- else\r
- rya = 0;\r
-\r
- if (yb[2] != 0)\r
- ryb = 2;\r
- else if (yb[1] != 0)\r
- ryb = 1;\r
- else\r
- ryb = 0;\r
-\r
- //kleinste Ecke finden und fuer \delta verwenden\r
-\r
- if (ya[rya] > 0) {\r
- if (yb[ryb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y0[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y3[i];\r
- }\r
- } else {\r
- if (yb[ryb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y1[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y2[i];\r
- //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
- }\r
- }\r
-\r
- // for (i = 0; i<3;++i)\r
- // d[i] = y0[i]-x0[i];\r
-\r
- //printf("(%d,%d)",rx,ry);\r
- if (rx == ry) {\r
- //printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
- //printf("[%.2f,%.2f@%.2f,%.2f]",xa[rxa],xb[rxb],ya[rya],yb[ryb]);\r
- //printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
- if (rxa == rya) {\r
- //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
- // y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
- tmp\r
- = calcParIntA( fabs(xa[rxa]), fabs(xb[rxb]),\r
- fabs(ya[rxa]), fabs(yb[rxb]), d[rxa],\r
- d[rxb], d[rx]);\r
- // printf("%d%d|%.2f\n",j,k,tmp);\r
- } else {\r
- //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
- // y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
- tmp\r
- = calcParIntA( fabs(xa[rxa]), fabs(xb[rxb]),\r
- fabs(yb[rxa]), fabs(ya[rxb]), d[rxa],\r
- d[rxb], d[rx]);\r
- // printf("%d%d|%.2f\n",j,k,tmp);\r
- }\r
-\r
- } else {\r
- //printf("(%d,%d@%d,%d)", rxa, rxb, rya, ryb);\r
- //printf("(%d,%d)", rx, ry);\r
- //printf("\n");\r
-\r
- if (rxa == rya) {\r
- tmp\r
- = calcOrtIntA( fabs(xb[rxb]), fabs(xa[rxa]),\r
- fabs(ya[rya]), fabs(yb[ryb]), d[rxb],\r
- d[rxa], d[rx]);\r
- } else if (rxa == ryb) {\r
- tmp\r
- = calcOrtIntA( fabs(xb[rxb]), fabs(xa[rxa]),\r
- fabs(yb[ryb]), fabs(ya[rya]), d[rxb],\r
- d[rxa], d[rx]);\r
- } else if (rxb == rya) {\r
- tmp\r
- = calcOrtIntA( fabs(xa[rxa]), fabs(xb[rxb]),\r
- fabs(ya[rya]), fabs(yb[ryb]), d[rxa],\r
- d[rxb], d[rx]);\r
- } else {\r
- tmp\r
- = calcOrtIntA( fabs(xa[rxa]), fabs(xb[rxb]),\r
- fabs(yb[ryb]), fabs(ya[rya]), d[rxa],\r
- d[rxb], d[rx]);\r
- }\r
-\r
- }\r
- A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
-\r
- }\r
-\r
- }\r
- //Rueckgabe (eventuell zurueck schreiben)\r
- // mxFree(x0);\r
- //mxFree(x1);\r
- //mxFree(x3);\r
- //mxFree(y0);\r
- //mxFree(y1);\r
- //mxFree(y3);\r
- //mxFree(xn);\r
- //mxFree(yn);\r
- //mxFree(d);\r
-\r
- return;\r
-}\r
--- /dev/null
+/*\r
+ * Art der Berechnung durch Parameter bestimmt...\r
+ * Analytisch/Semianalytisch/...\r
+ *\r
+ *\r
+ */\r
+\r
+#include <iostream>\r
+#include <cmath>\r
+#include <cassert>\r
+#include "mex.h"\r
+#include <stdlib.h>\r
+\r
+#define M_EPS 10e-8\r
+\r
+//#include "tbb/parallel_for.h"\r
+\r
+\r
+#include "slpRectangle.hpp"\r
+\r
+using namespace std;\r
+\r
+int dimOfVec(double* vec) {\r
+ if (fabs(vec[2]) > M_EPS)\r
+ return 2;\r
+ else if (fabs(vec[1]) > M_EPS)\r
+ return 1;\r
+ else if (fabs(vec[0]) > M_EPS)\r
+ return 0;\r
+ else {\r
+ cerr << "dimOfVec : (" << vec[0] << vec[1] << vec[2]\r
+ << ") alle Werte 0 \n";\r
+ return -1;\r
+ }\r
+\r
+}\r
+\r
+inline int dimOfThird(int a, int b) {\r
+ return ((-(a + b) % 3) + 3) % 3;\r
+}\r
+\r
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
+\r
+ int i, j, k; //Schleifenindizes\r
+ double tmp; //Zwischenspeicherung der Einzelnen Werte\r
+\r
+ //Sicherheitsabfragen zu Datengroessen\r
+ if ((nrhs != 4))\r
+ mexErrMsgTxt(\r
+ "expected (coordinates(Nx3),elements(Mx4),mu(double),type(int))");\r
+ if (nlhs > 1)\r
+ mexErrMsgTxt("has only one output argument");\r
+\r
+ int cm = mxGetM(prhs[0]);\r
+ int cn = mxGetN(prhs[0]);\r
+ if (cn != 3)\r
+ mexErrMsgTxt("expected coordinates (Nx3)");\r
+ cout << " Coordinaten:" << cm << endl;\r
+\r
+ int em = mxGetM(prhs[1]);\r
+ int en = mxGetN(prhs[1]);\r
+ if (en != 4)\r
+ mexErrMsgTxt("expected elements (Mx4)");\r
+ cout << " Elemente:" << em << endl;\r
+\r
+ //Auslesen der Parameter\r
+\r
+ plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);\r
+ double * A = mxGetPr(plhs[0]);\r
+ double * C = mxGetPr(prhs[0]);\r
+ double * E = mxGetPr(prhs[1]);\r
+\r
+ double mu = *(mxGetPr(prhs[2]));\r
+ int type = *(mxGetPr(prhs[3]));\r
+\r
+ //Initialisieren der Hilfsvektoren\r
+\r
+ double * x0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+\r
+ double * y0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+\r
+ double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+ double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
+\r
+ // Welche Funktion soll verwendet werden\r
+ double(*ctypeP)(double, double, double, double, double, double, double,\r
+ double);\r
+ double(*ctypeO)(double, double, double, double, double, double, double,\r
+ double);\r
+\r
+ //Art der Berechnung bestimmen\r
+ cout << " Type:" << type << endl;\r
+ switch (type) {\r
+ case 0:\r
+ ctypeP = calcParIntO0;\r
+ ctypeO = calcOrtIntO0;\r
+ case 1:\r
+ ctypeP = calcParIntO1;\r
+ ctypeO = calcOrtIntO1;\r
+ case 2:\r
+ ctypeP = calcParIntO2;\r
+ ctypeO = calcOrtIntO2;\r
+ case 3:\r
+ ctypeP = calcParIntO3;\r
+ ctypeO = calcOrtIntO3;\r
+\r
+ }\r
+\r
+ //LageInformationen\r
+ int rx, rxa, rxb, ry, rya, ryb;\r
+\r
+ //Ausrechnen\r
+ for (j = 0; j < em; ++j) {\r
+ x0[0] = C[(int) E[j] - 1];\r
+ x0[1] = C[cm + (int) E[j] - 1];\r
+ x0[2] = C[2 * cm + (int) E[j] - 1];\r
+\r
+ x1[0] = C[(int) E[em + j] - 1];\r
+ x1[1] = C[cm + (int) E[em + j] - 1];\r
+ x1[2] = C[2 * cm + (int) E[em + j] - 1];\r
+\r
+ x2[0] = C[(int) E[2 * em + j] - 1];\r
+ x2[1] = C[cm + (int) E[2 * em + j] - 1];\r
+ x2[2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
+\r
+ x3[0] = C[(int) E[3 * em + j] - 1];\r
+ x3[1] = C[cm + (int) E[3 * em + j] - 1];\r
+ x3[2] = C[2 * cm + (int) E[3 * em + j] - 1];\r
+\r
+ for (i = 0; i < 3; ++i)\r
+ xa[i] = x1[i] - x0[i];\r
+\r
+ for (i = 0; i < 3; ++i)\r
+ xb[i] = x3[i] - x0[i];\r
+\r
+ // Lageeigenschaften des Flaechenstuecks\r
+\r
+ rxa = dimOfVec(xa);\r
+ rxb = dimOfVec(xb);\r
+ rx = dimOfThird(rxa, rxb);\r
+\r
+ //kleinste Ecke finden und fuer \delta verwenden\r
+\r
+ if (xa[rxa] > 0) {\r
+ if (xb[rxb] > 0) {\r
+ for (i = 0; i < 3; ++i)\r
+ dt[i] = -x0[i];\r
+ } else {\r
+ for (i = 0; i < 3; ++i)\r
+ dt[i] = -x3[i];\r
+ }\r
+ } else {\r
+ if (xb[rxb] > 0) {\r
+ for (i = 0; i < 3; ++i)\r
+ dt[i] = -x1[i];\r
+ } else {\r
+ for (i = 0; i < 3; ++i)\r
+ dt[i] = -x2[i];\r
+ }\r
+ }\r
+\r
+ for (k = 0; k < em; ++k) {\r
+ y0[0] = C[(int) E[k] - 1];\r
+ y0[1] = C[cm + (int) E[k] - 1];\r
+ y0[2] = C[2 * cm + (int) E[k] - 1];\r
+\r
+ y1[0] = C[(int) E[em + k] - 1];\r
+ y1[1] = C[cm + (int) E[em + k] - 1];\r
+ y1[2] = C[2 * cm + (int) E[em + k] - 1];\r
+\r
+ y2[0] = C[(int) E[2 * em + k] - 1];\r
+ y2[1] = C[cm + (int) E[2 * em + k] - 1];\r
+ y2[2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
+\r
+ y3[0] = C[(int) E[3 * em + k] - 1];\r
+ y3[1] = C[cm + (int) E[3 * em + k] - 1];\r
+ y3[2] = C[2 * cm + (int) E[3 * em + k] - 1];\r
+\r
+ for (i = 0; i < 3; ++i)\r
+ ya[i] = y1[i] - y0[i];\r
+\r
+ for (i = 0; i < 3; ++i)\r
+ yb[i] = y3[i] - y0[i];\r
+\r
+ // Lageeigenschaften des Flaechenstuecks\r
+\r
+ rya = dimOfVec(ya);\r
+ ryb = dimOfVec(yb);\r
+ ry = dimOfThird(rya, ryb);\r
+\r
+ //kleinste Ecke finden und fuer \delta verwenden\r
+\r
+ if (ya[rya] > 0) {\r
+ if (yb[ryb] > 0) {\r
+ for (i = 0; i < 3; ++i)\r
+ d[i] = dt[i] + y0[i];\r
+ } else {\r
+ for (i = 0; i < 3; ++i)\r
+ d[i] = dt[i] + y3[i];\r
+ }\r
+ } else {\r
+ if (yb[ryb] > 0) {\r
+ for (i = 0; i < 3; ++i)\r
+ d[i] = dt[i] + y1[i];\r
+ } else {\r
+ for (i = 0; i < 3; ++i)\r
+ d[i] = dt[i] + y2[i];\r
+ }\r
+ }\r
+\r
+ if (rx == ry) {\r
+ if (rxa == rya) {\r
+ tmp\r
+ = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]),\r
+ fabs(ya[rxa]), fabs(yb[rxb]), d[rxa],\r
+ d[rxb], d[rx], mu);\r
+\r
+ } else {\r
+ tmp\r
+ = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]),\r
+ fabs(yb[rxa]), fabs(ya[rxb]), d[rxa],\r
+ d[rxb], d[rx], mu);\r
+ }\r
+\r
+ } else {\r
+\r
+ if (rxa == rya) {\r
+ tmp\r
+ = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]),\r
+ fabs(ya[rya]), fabs(yb[ryb]), d[rxb],\r
+ d[rxa], d[rx], mu);\r
+ } else if (rxa == ryb) {\r
+ tmp\r
+ = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]),\r
+ fabs(yb[ryb]), fabs(ya[rya]), d[rxb],\r
+ d[rxa], d[rx], mu);\r
+ } else if (rxb == rya) {\r
+ tmp\r
+ = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]),\r
+ fabs(ya[rya]), fabs(yb[ryb]), d[rxa],\r
+ d[rxb], d[rx], mu);\r
+ } else {\r
+ tmp\r
+ = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]),\r
+ fabs(yb[ryb]), fabs(ya[rya]), d[rxa],\r
+ d[rxb], d[rx], mu);\r
+ }\r
+\r
+ }\r
+ A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
+\r
+ }\r
+\r
+ }\r
+ //Rueckgabe (eventuell zurueck schreiben)\r
+ //mxFree(x0);\r
+ //mxFree(x1);\r
+ //mxFree(x3);\r
+ //mxFree(y0);\r
+ //mxFree(y1);\r
+ //mxFree(y3);\r
+ //mxFree(d);\r
+ //mxFree(dt);\r
+\r
+ return;\r
+}\r
+++ /dev/null
-/*\r
- * voll Analytisch\r
- * Element mit groesserer Diagonale nach aussen gedreht\r
- *\r
- *\r
- */\r
-\r
-\r
-#include <iostream>\r
-#include <cmath>\r
-#include <cassert>\r
-#include "mex.h"\r
-#include <stdlib.h>\r
-\r
-//#include "tbb/parallel_for.h"\r
-\r
-\r
-#include "slpRectangle.hpp"\r
-\r
-#define EPS 0.02\r
-\r
-using namespace std;\r
-\r
-inline int dimOfVec(double * vec)\r
-{\r
- if (vec[2] != 0)\r
- return 2;\r
- else if (vec[1] != 0)\r
- return 1;\r
- else\r
- return 0;\r
-}\r
-\r
-\r
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
-\r
- int i, j, k;\r
-\r
- //sicherheitsabfragen zu Datengroessen\r
- if (nrhs != 2)\r
- mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx4))");\r
- if (nlhs > 1)\r
- mexErrMsgTxt("has only one output argument");\r
-\r
- int cm = mxGetM(prhs[0]);\r
- int cn = mxGetN(prhs[0]);\r
-\r
- if (cn != 3)\r
- mexErrMsgTxt("expected coordinates (Nx3)");\r
- int em = mxGetM(prhs[1]);\r
- int en = mxGetN(prhs[1]);\r
- if (en != 4)\r
- mexErrMsgTxt("expected elements (Mx4)");\r
- //Vorbereitung der Daten\r
-\r
- plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);\r
- double * A = mxGetPr(plhs[0]);\r
- double * C = mxGetPr(prhs[0]);\r
- double * E = mxGetPr(prhs[1]);\r
-\r
- double * x0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double * y0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * yn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double tmp,tmp2;\r
-\r
- int itmp;\r
- double dtmp;\r
- double * vtmp;\r
-\r
- long double xda,xdb,yda,ydb;\r
-\r
- //LageInformationen\r
- int rx, rxa, rxb, ry, rya, ryb;\r
-\r
- //Ausrechnen\r
- for (j = 0; j < em; ++j) {\r
- x0[0] = C[(int) E[j] - 1];\r
- x0[1] = C[cm + (int) E[j] - 1];\r
- x0[2] = C[2 * cm + (int) E[j] - 1];\r
-\r
- x1[0] = C[(int) E[em + j] - 1];\r
- x1[1] = C[cm + (int) E[em + j] - 1];\r
- x1[2] = C[2 * cm + (int) E[em + j] - 1];\r
-\r
- x2[0] = C[(int) E[2 * em + j] - 1];\r
- x2[1] = C[cm + (int) E[2 * em + j] - 1];\r
- x2[2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
-\r
- x3[0] = C[(int) E[3 * em + j] - 1];\r
- x3[1] = C[cm + (int) E[3 * em + j] - 1];\r
- x3[2] = C[2 * cm + (int) E[3 * em + j] - 1];\r
-\r
- for (i = 0; i < 3; ++i)\r
- xa[i] = x1[i] - x0[i];\r
-\r
- for (i = 0; i < 3; ++i)\r
- xb[i] = x3[i] - x0[i];\r
-\r
- //Kreuzprodukt\r
- xn[0] = (xa[1]) * (xb[2]) - (xa[2]) * (xb[1]);\r
- xn[1] = (xa[2]) * (xb[0]) - (xa[0]) * (xb[2]);\r
- xn[2] = (xa[0]) * (xb[1]) - (xa[1]) * (xb[0]);\r
-\r
- // Lageeigenschaften des Flaechenstuecks\r
- //if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])\r
-\r
- rx = dimOfVec(xn);\r
-\r
- rxa = dimOfVec(xa);\r
-\r
- rxb = dimOfVec(xb);\r
-\r
-\r
- //kleinste Ecke finden und fuer \delta verwenden\r
-\r
- if (xa[rxa] > 0) {\r
- if (xb[rxb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x0[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x3[i];\r
- }\r
- } else {\r
- if (xb[rxb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x1[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x2[i];\r
- //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
- }\r
- }\r
-\r
- xda = xa[rxa];\r
- xdb = xb[rxb];\r
-\r
- for (k = 0; k < em; ++k) {\r
-\r
-\r
-\r
- y0[0] = C[(int) E[k] - 1];\r
- y0[1] = C[cm + (int) E[k] - 1];\r
- y0[2] = C[2 * cm + (int) E[k] - 1];\r
-\r
- y1[0] = C[(int) E[em + k] - 1];\r
- y1[1] = C[cm + (int) E[em + k] - 1];\r
- y1[2] = C[2 * cm + (int) E[em + k] - 1];\r
-\r
- y2[0] = C[(int) E[2 * em + k] - 1];\r
- y2[1] = C[cm + (int) E[2 * em + k] - 1];\r
- y2[2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
-\r
- y3[0] = C[(int) E[3 * em + k] - 1];\r
- y3[1] = C[cm + (int) E[3 * em + k] - 1];\r
- y3[2] = C[2 * cm + (int) E[3 * em + k] - 1];\r
-\r
- for (i = 0; i < 3; ++i)\r
- ya[i] = y1[i] - y0[i];\r
-\r
- for (i = 0; i < 3; ++i)\r
- yb[i] = y3[i] - y0[i];\r
-\r
- yn[0] = (ya[1]) * (yb[2]) - (ya[2]) * (yb[1]);\r
- yn[1] = (ya[2]) * (yb[0]) - (ya[0]) * (yb[2]);\r
- yn[2] = (ya[0]) * (yb[1]) - (ya[1]) * (yb[0]);\r
-\r
- ry = dimOfVec(yn);\r
-\r
- rya = dimOfVec(ya);\r
-\r
- ryb = dimOfVec(yb);\r
-\r
- //kleinste Ecke finden und fuer \delta verwenden\r
-\r
- if (ya[rya] > 0) {\r
- if (yb[ryb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y0[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y3[i];\r
- }\r
- } else {\r
- if (yb[ryb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y1[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y2[i];\r
- //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
- }\r
- }\r
-\r
- yda = ya[rya];\r
- ydb = yb[ryb];\r
-\r
- //printf("(%d,%d)",rx,ry);\r
- if (rx == ry) {\r
- //printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
- //printf("[%.2f,%.2f@%.2f,%.2f]",xda,xdb,yda,ydb);\r
- //printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
- if (rxa == rya) {\r
- //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
- // y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
-\r
- tmp\r
- = calcParIntO1( fabs(xda), fabs(xdb),\r
- fabs(yda), fabs(ydb), d[rxa],\r
- d[rxb], d[rx]);\r
-\r
- } else {\r
- //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
- // y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
- tmp\r
- = calcParIntO1( fabs(xda), fabs(xdb),\r
- fabs(ydb), fabs(yda), d[rxa],\r
- d[rxb], d[rx]);\r
- // printf("%d%d|%.2f\n",j,k,tmp);\r
- }\r
-\r
- } else {\r
- //printf("(%d,%d@%d,%d)", rxa, rxb, rya, ryb);\r
- //printf("(%d,%d)", rx, ry);\r
- //printf("\n");\r
-\r
- if (rxa == rya) {\r
- tmp\r
- = calcOrtIntO1( fabs(xdb), fabs(xda),\r
- fabs(yda), fabs(ydb), d[rxb],\r
- d[rxa], d[rx]);\r
- } else if (rxa == ryb) {\r
- tmp\r
- = calcOrtIntO1( fabs(xdb), fabs(xda),\r
- fabs(ydb), fabs(yda), d[rxb],\r
- d[rxa], d[rx]);\r
- } else if (rxb == rya) {\r
- tmp\r
- = calcOrtIntO1( fabs(xda), fabs(xdb),\r
- fabs(yda), fabs(ydb), d[rxa],\r
- d[rxb], d[rx]);\r
- } else {\r
- tmp\r
- = calcOrtIntO1( fabs(xda), fabs(xdb),\r
- fabs(ydb), fabs(yda), d[rxa],\r
- d[rxb], d[rx]);\r
- }\r
-\r
- }\r
- A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
-\r
- }\r
-\r
- }\r
- //printf("\n");\r
- //Rueckgabe (eventuell zurueck schreiben)\r
- // mxFree(x0);\r
- //mxFree(x1);\r
- //mxFree(x3);\r
- //mxFree(y0);\r
- //mxFree(y1);\r
- //mxFree(y3);\r
- //mxFree(xn);\r
- //mxFree(yn);\r
- //mxFree(d);\r
-\r
- return;\r
-}\r
+++ /dev/null
-/*\r
- * voll Analytisch\r
- * Element mit groesserer Diagonale nach aussen gedreht\r
- * Quadratur ueber E_j bei dist(E_j,E_k) >= eta*dia(E_j) | dia(E_j) <= dia(E_k)\r
- *\r
- */\r
-\r
-\r
-#include <iostream>\r
-#include <cmath>\r
-#include <cassert>\r
-#include "mex.h"\r
-#include <stdlib.h>\r
-\r
-//#include "tbb/parallel_for.h"\r
-\r
-\r
-#include "slpRectangle.hpp"\r
-\r
-//Quadraturverwendung!!!\r
-#define eta 0.8\r
-\r
-using namespace std;\r
-\r
-inline int dimOfVec(double * vec)\r
-{\r
-\r
- if (vec[2] != 0)\r
- return 2;\r
- else if (vec[1] != 0)\r
- return 1;\r
- else if (vec[0] != 0)\r
- return 0;\r
- else{\r
- cerr << "dimOfVec : alle Werte 0";\r
- return -1;\r
- }\r
-}\r
-\r
-inline int dimOfThird(int a, int b)\r
-{\r
- return mod(-(a+b),3);\r
-}\r
-\r
-\r
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
-\r
- int i, j, k;\r
-\r
- //sicherheitsabfragen zu Datengroessen\r
- if (nrhs != 2)\r
- mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx4))");\r
- if (nlhs > 1)\r
- mexErrMsgTxt("has only one output argument");\r
-\r
- int cm = mxGetM(prhs[0]);\r
- int cn = mxGetN(prhs[0]);\r
-\r
- if (cn != 3)\r
- mexErrMsgTxt("expected coordinates (Nx3)");\r
- int em = mxGetM(prhs[1]);\r
- int en = mxGetN(prhs[1]);\r
- if (en != 4)\r
- mexErrMsgTxt("expected elements (Mx4)");\r
- //Vorbereitung der Daten\r
-\r
- plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);\r
- double * A = mxGetPr(plhs[0]);\r
- double * C = mxGetPr(prhs[0]);\r
- double * E = mxGetPr(prhs[1]);\r
-\r
- double * x0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double * y0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * yn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double tmp,tmp2;\r
-\r
- int itmp;\r
- double dtmp;\r
- double * vtmp;\r
-\r
- long double xda,xdb,yda,ydb;\r
-\r
- //LageInformationen\r
- int rx, rxa, rxb, ry, rya, ryb;\r
-\r
- //Ausrechnen\r
- for (j = 0; j < em; ++j) {\r
- x0[0] = C[(int) E[j] - 1];\r
- x0[1] = C[cm + (int) E[j] - 1];\r
- x0[2] = C[2 * cm + (int) E[j] - 1];\r
-\r
- x1[0] = C[(int) E[em + j] - 1];\r
- x1[1] = C[cm + (int) E[em + j] - 1];\r
- x1[2] = C[2 * cm + (int) E[em + j] - 1];\r
-\r
- x2[0] = C[(int) E[2 * em + j] - 1];\r
- x2[1] = C[cm + (int) E[2 * em + j] - 1];\r
- x2[2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
-\r
- x3[0] = C[(int) E[3 * em + j] - 1];\r
- x3[1] = C[cm + (int) E[3 * em + j] - 1];\r
- x3[2] = C[2 * cm + (int) E[3 * em + j] - 1];\r
-\r
- for (i = 0; i < 3; ++i)\r
- xa[i] = x1[i] - x0[i];\r
-\r
- for (i = 0; i < 3; ++i)\r
- xb[i] = x3[i] - x0[i];\r
-\r
- //Kreuzprodukt\r
- xn[0] = (xa[1]) * (xb[2]) - (xa[2]) * (xb[1]);\r
- xn[1] = (xa[2]) * (xb[0]) - (xa[0]) * (xb[2]);\r
- xn[2] = (xa[0]) * (xb[1]) - (xa[1]) * (xb[0]);\r
-\r
- // Lageeigenschaften des Flaechenstuecks\r
- //if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])\r
-\r
- rx = dimOfVec(xn);\r
-\r
- rxa = dimOfVec(xa);\r
-\r
- rxb = dimOfVec(xb);\r
-\r
-\r
- //kleinste Ecke finden und fuer \delta verwenden\r
-\r
- if (xa[rxa] > 0) {\r
- if (xb[rxb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x0[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x3[i];\r
- }\r
- } else {\r
- if (xb[rxb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x1[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x2[i];\r
- //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
- }\r
- }\r
-\r
- xda = xa[rxa];\r
- xdb = xb[rxb];\r
-\r
- for (k = 0; k < em; ++k) {\r
-\r
-\r
-\r
- y0[0] = C[(int) E[k] - 1];\r
- y0[1] = C[cm + (int) E[k] - 1];\r
- y0[2] = C[2 * cm + (int) E[k] - 1];\r
-\r
- y1[0] = C[(int) E[em + k] - 1];\r
- y1[1] = C[cm + (int) E[em + k] - 1];\r
- y1[2] = C[2 * cm + (int) E[em + k] - 1];\r
-\r
- y2[0] = C[(int) E[2 * em + k] - 1];\r
- y2[1] = C[cm + (int) E[2 * em + k] - 1];\r
- y2[2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
-\r
- y3[0] = C[(int) E[3 * em + k] - 1];\r
- y3[1] = C[cm + (int) E[3 * em + k] - 1];\r
- y3[2] = C[2 * cm + (int) E[3 * em + k] - 1];\r
-\r
- for (i = 0; i < 3; ++i)\r
- ya[i] = y1[i] - y0[i];\r
-\r
- for (i = 0; i < 3; ++i)\r
- yb[i] = y3[i] - y0[i];\r
-\r
- yn[0] = (ya[1]) * (yb[2]) - (ya[2]) * (yb[1]);\r
- yn[1] = (ya[2]) * (yb[0]) - (ya[0]) * (yb[2]);\r
- yn[2] = (ya[0]) * (yb[1]) - (ya[1]) * (yb[0]);\r
-\r
- ry = dimOfVec(yn);\r
-\r
- rya = dimOfVec(ya);\r
-\r
- ryb = dimOfVec(yb);\r
-\r
- //kleinste Ecke finden und fuer \delta verwenden\r
-\r
- if (ya[rya] > 0) {\r
- if (yb[ryb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y0[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y3[i];\r
- }\r
- } else {\r
- if (yb[ryb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y1[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y2[i];\r
- //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
- }\r
- }\r
-\r
- yda = ya[rya];\r
- ydb = yb[ryb];\r
-\r
- //printf("(%d,%d)",rx,ry);\r
- if (rx == ry) {\r
- //printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
- //printf("[%.2f,%.2f@%.2f,%.2f]",xda,xdb,yda,ydb);\r
- //printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
- if (rxa == rya) {\r
- //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
- // y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
-\r
- tmp\r
- = calcParIntO2( fabs(xda), fabs(xdb),\r
- fabs(yda), fabs(ydb), d[rxa],\r
- d[rxb], d[rx],eta);\r
-\r
- } else {\r
- //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
- // y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
- tmp\r
- = calcParIntO2( fabs(xda), fabs(xdb),\r
- fabs(ydb), fabs(yda), d[rxa],\r
- d[rxb], d[rx],eta);\r
- // printf("%d%d|%.2f\n",j,k,tmp);\r
- }\r
-\r
- } else {\r
- //printf("(%d,%d@%d,%d)", rxa, rxb, rya, ryb);\r
- //printf("(%d,%d)", rx, ry);\r
- //printf("\n");\r
-\r
- if (rxa == rya) {\r
- tmp\r
- = calcOrtIntO2( fabs(xdb), fabs(xda),\r
- fabs(yda), fabs(ydb), d[rxb],\r
- d[rxa], d[rx],eta);\r
- } else if (rxa == ryb) {\r
- tmp\r
- = calcOrtIntO2( fabs(xdb), fabs(xda),\r
- fabs(ydb), fabs(yda), d[rxb],\r
- d[rxa], d[rx],eta);\r
- } else if (rxb == rya) {\r
- tmp\r
- = calcOrtIntO2( fabs(xda), fabs(xdb),\r
- fabs(yda), fabs(ydb), d[rxa],\r
- d[rxb], d[rx],eta);\r
- } else {\r
- tmp\r
- = calcOrtIntO2( fabs(xda), fabs(xdb),\r
- fabs(ydb), fabs(yda), d[rxa],\r
- d[rxb], d[rx],eta);\r
- }\r
-\r
- }\r
- A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
-\r
- }\r
-\r
- }\r
- //printf("\n");\r
- //Rueckgabe (eventuell zurueck schreiben)\r
- // mxFree(x0);\r
- //mxFree(x1);\r
- //mxFree(x3);\r
- //mxFree(y0);\r
- //mxFree(y1);\r
- //mxFree(y3);\r
- //mxFree(xn);\r
- //mxFree(yn);\r
- //mxFree(d);\r
-\r
- return;\r
-}\r
+++ /dev/null
-/*\r
- * voll Analytisch\r
- * Element mit groesserer Diagonale nach aussen gedreht\r
- * Quadratur ueber E_j bei dist(E_j,E_k) >= eta*dia(E_j) | dia(E_j) <= dia(E_k)\r
- *\r
- */\r
-\r
-\r
-#include <iostream>\r
-#include <cmath>\r
-#include <cassert>\r
-#include "mex.h"\r
-#include <stdlib.h>\r
-\r
-//#include "tbb/parallel_for.h"\r
-\r
-\r
-#include "slpRectangle.hpp"\r
-\r
-//Quadraturverwendung!!!\r
-#define eta 0.8\r
-\r
-using namespace std;\r
-\r
-inline int dimOfVec(double * vec)\r
-{\r
- if (vec[2] != 0)\r
- return 2;\r
- else if (vec[1] != 0)\r
- return 1;\r
- else\r
- return 0;\r
-}\r
-\r
-\r
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
-\r
- int i, j, k;\r
-\r
- //sicherheitsabfragen zu Datengroessen\r
- if (nrhs != 2)\r
- mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx4))");\r
- if (nlhs > 1)\r
- mexErrMsgTxt("has only one output argument");\r
-\r
- int cm = mxGetM(prhs[0]);\r
- int cn = mxGetN(prhs[0]);\r
-\r
- if (cn != 3)\r
- mexErrMsgTxt("expected coordinates (Nx3)");\r
- int em = mxGetM(prhs[1]);\r
- int en = mxGetN(prhs[1]);\r
- if (en != 4)\r
- mexErrMsgTxt("expected elements (Mx4)");\r
- //Vorbereitung der Daten\r
-\r
- plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);\r
- double * A = mxGetPr(plhs[0]);\r
- double * C = mxGetPr(prhs[0]);\r
- double * E = mxGetPr(prhs[1]);\r
-\r
- double * x0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double * y0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * yn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
- double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double tmp,tmp2;\r
-\r
- int itmp;\r
- double dtmp;\r
- double * vtmp;\r
-\r
- long double xda,xdb,yda,ydb;\r
-\r
- //LageInformationen\r
- int rx, rxa, rxb, ry, rya, ryb;\r
-\r
- //Ausrechnen\r
- for (j = 0; j < em; ++j) {\r
- x0[0] = C[(int) E[j] - 1];\r
- x0[1] = C[cm + (int) E[j] - 1];\r
- x0[2] = C[2 * cm + (int) E[j] - 1];\r
-\r
- x1[0] = C[(int) E[em + j] - 1];\r
- x1[1] = C[cm + (int) E[em + j] - 1];\r
- x1[2] = C[2 * cm + (int) E[em + j] - 1];\r
-\r
- x2[0] = C[(int) E[2 * em + j] - 1];\r
- x2[1] = C[cm + (int) E[2 * em + j] - 1];\r
- x2[2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
-\r
- x3[0] = C[(int) E[3 * em + j] - 1];\r
- x3[1] = C[cm + (int) E[3 * em + j] - 1];\r
- x3[2] = C[2 * cm + (int) E[3 * em + j] - 1];\r
-\r
- for (i = 0; i < 3; ++i)\r
- xa[i] = x1[i] - x0[i];\r
-\r
- for (i = 0; i < 3; ++i)\r
- xb[i] = x3[i] - x0[i];\r
-\r
- //Kreuzprodukt\r
- xn[0] = (xa[1]) * (xb[2]) - (xa[2]) * (xb[1]);\r
- xn[1] = (xa[2]) * (xb[0]) - (xa[0]) * (xb[2]);\r
- xn[2] = (xa[0]) * (xb[1]) - (xa[1]) * (xb[0]);\r
-\r
- // Lageeigenschaften des Flaechenstuecks\r
- //if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])\r
-\r
- rx = dimOfVec(xn);\r
-\r
- rxa = dimOfVec(xa);\r
-\r
- rxb = dimOfVec(xb);\r
-\r
-\r
- //kleinste Ecke finden und fuer \delta verwenden\r
-\r
- if (xa[rxa] > 0) {\r
- if (xb[rxb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x0[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x3[i];\r
- }\r
- } else {\r
- if (xb[rxb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x1[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- dt[i] = -x2[i];\r
- //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
- }\r
- }\r
-\r
- xda = xa[rxa];\r
- xdb = xb[rxb];\r
-\r
- for (k = 0; k < em; ++k) {\r
-\r
-\r
-\r
- y0[0] = C[(int) E[k] - 1];\r
- y0[1] = C[cm + (int) E[k] - 1];\r
- y0[2] = C[2 * cm + (int) E[k] - 1];\r
-\r
- y1[0] = C[(int) E[em + k] - 1];\r
- y1[1] = C[cm + (int) E[em + k] - 1];\r
- y1[2] = C[2 * cm + (int) E[em + k] - 1];\r
-\r
- y2[0] = C[(int) E[2 * em + k] - 1];\r
- y2[1] = C[cm + (int) E[2 * em + k] - 1];\r
- y2[2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
-\r
- y3[0] = C[(int) E[3 * em + k] - 1];\r
- y3[1] = C[cm + (int) E[3 * em + k] - 1];\r
- y3[2] = C[2 * cm + (int) E[3 * em + k] - 1];\r
-\r
- for (i = 0; i < 3; ++i)\r
- ya[i] = y1[i] - y0[i];\r
-\r
- for (i = 0; i < 3; ++i)\r
- yb[i] = y3[i] - y0[i];\r
-\r
- yn[0] = (ya[1]) * (yb[2]) - (ya[2]) * (yb[1]);\r
- yn[1] = (ya[2]) * (yb[0]) - (ya[0]) * (yb[2]);\r
- yn[2] = (ya[0]) * (yb[1]) - (ya[1]) * (yb[0]);\r
-\r
- ry = dimOfVec(yn);\r
-\r
- rya = dimOfVec(ya);\r
-\r
- ryb = dimOfVec(yb);\r
-\r
- //kleinste Ecke finden und fuer \delta verwenden\r
-\r
- if (ya[rya] > 0) {\r
- if (yb[ryb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y0[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y3[i];\r
- }\r
- } else {\r
- if (yb[ryb] > 0) {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y1[i];\r
- } else {\r
- for (i = 0; i < 3; ++i)\r
- d[i] = dt[i] + y2[i];\r
- //printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");\r
- }\r
- }\r
-\r
- yda = ya[rya];\r
- ydb = yb[ryb];\r
-\r
- //printf("(%d,%d)",rx,ry);\r
- if (rx == ry) {\r
- //printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);\r
- //printf("[%.2f,%.2f@%.2f,%.2f]",xda,xdb,yda,ydb);\r
- //printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
- if (rxa == rya) {\r
- //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
- // y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
-\r
- tmp\r
- = calcParIntO3( fabs(xda), fabs(xdb),\r
- fabs(yda), fabs(ydb), d[rxa],\r
- d[rxb], d[rx],eta);\r
-\r
- } else {\r
- //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],\r
- // y1[0], y0[1], y3[1], d[0], d[1], d[2]);\r
- //printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);\r
- tmp\r
- = calcParIntO3( fabs(xda), fabs(xdb),\r
- fabs(ydb), fabs(yda), d[rxa],\r
- d[rxb], d[rx],eta);\r
- // printf("%d%d|%.2f\n",j,k,tmp);\r
- }\r
-\r
- } else {\r
- //printf("(%d,%d@%d,%d)", rxa, rxb, rya, ryb);\r
- //printf("(%d,%d)", rx, ry);\r
- //printf("\n");\r
-\r
- if (rxa == rya) {\r
- tmp\r
- = calcOrtIntO3( fabs(xdb), fabs(xda),\r
- fabs(yda), fabs(ydb), d[rxb],\r
- d[rxa], d[rx],eta);\r
- } else if (rxa == ryb) {\r
- tmp\r
- = calcOrtIntO3( fabs(xdb), fabs(xda),\r
- fabs(ydb), fabs(yda), d[rxb],\r
- d[rxa], d[rx],eta);\r
- } else if (rxb == rya) {\r
- tmp\r
- = calcOrtIntO3( fabs(xda), fabs(xdb),\r
- fabs(yda), fabs(ydb), d[rxa],\r
- d[rxb], d[rx],eta);\r
- } else {\r
- tmp\r
- = calcOrtIntO3( fabs(xda), fabs(xdb),\r
- fabs(ydb), fabs(yda), d[rxa],\r
- d[rxb], d[rx],eta);\r
- }\r
-\r
- }\r
- A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
-\r
- }\r
-\r
- }\r
- //printf("\n");\r
- //Rueckgabe (eventuell zurueck schreiben)\r
- // mxFree(x0);\r
- //mxFree(x1);\r
- //mxFree(x3);\r
- //mxFree(y0);\r
- //mxFree(y1);\r
- //mxFree(y3);\r
- //mxFree(xn);\r
- //mxFree(yn);\r
- //mxFree(d);\r
-\r
- return;\r
-}\r
+++ /dev/null
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.hpp"
-
-#define my_PI 3.141592653589793
-#define EPS 0.02
-
-using namespace std;
-
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
- //sicherheitsabfragen zu Datengroessen
- if (nrhs != 7)
- mexErrMsgTxt("expected (x1,x2,y1,y2,d1,d2,d3)");
- if (nlhs > 2)
- mexErrMsgTxt("has maximal two output arguments");
-
- plhs[0] = mxCreateDoubleMatrix(1, 3, mxREAL);
-
- *mxGetPr(plhs[0]) = calcortIntA(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-// *(mxGetPr(plhs[0])+1) = calcParIntQX1X2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-// *(mxGetPr(plhs[0])+2) = calcParIntQY1X1(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
-}
-
-
-
+++ /dev/null
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.hpp"
-
-#define my_PI 3.141592653589793
-#define EPS 0.02
-
-using namespace std;
-
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
- //sicherheitsabfragen zu Datengroessen
- if (nrhs != 7)
- mexErrMsgTxt("expected (x1,x2,y1,y2,d1,d2,d3)");
- if (nlhs > 2)
- mexErrMsgTxt("has maximal two output arguments");
-
-
- plhs[0] = mxCreateDoubleMatrix(1, 3, mxREAL);
-
- *mxGetPr(plhs[0]) = calcParIntA(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
- *(mxGetPr(plhs[0])+1) = calcParIntQX1X2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
- *(mxGetPr(plhs[0])+2) = calcParIntQY1X1(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
-
-}
-
-
-
+++ /dev/null
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.hpp"
-
-#define my_PI 3.141592653589793
-#define EPS 0.02
-
-using namespace std;
-
-
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
- //sicherheitsabfragen zu Datengroessen
- if (nrhs != 4)
- mexErrMsgTxt("expected (p, y, x, l)");
- if (nlhs > 2)
- mexErrMsgTxt("has maximal two output arguments");
-
-
- plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
- plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-
- *mxGetPr(plhs[0]) = g_AY(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
- *mxGetPr(plhs[0]) = g_QY(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
-
-}
-
-
-
+++ /dev/null
-
-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 mex_build_A.cpp slpRectangle.cpp
-mex mex_build_As1.cpp slpRectangle.cpp
-mex mex_build_As2.cpp slpRectangle.cpp
-
\r
#include "slpRectangle.hpp"\r
\r
-#define EPS 10e5\r
-\r
using namespace std;\r
\r
-double gauss_w5[] = {0.1185,0.2393,0.2844,0.2393,0.1185};\r
-double gauss_n5[] = {0.0469,0.2308,0.5000,0.7692,0.9531};\r
-\r
-\r
+double gauss_w5[] = { 0.1185, 0.2393, 0.2844, 0.2393, 0.1185 };\r
+double gauss_n5[] = { 0.0469, 0.2308, 0.5000, 0.7692, 0.9531 };\r
\r
int sign(double);\r
//double arsinh(double);\r
\r
/* ============================== */\r
\r
+\r
int inline sign(double x) {\r
return x > 0 ? 1 : (x < 0 ? -1 : 0);\r
}\r
}\r
*/\r
\r
-double inline max(double x, double y){\r
- return x<y?y:x;\r
+double inline max(double x, double y) {\r
+ return x < y ? y : x;\r
}\r
\r
double g_QY(double p, double y, double x, double l) {\r
\r
double sol = 0;\r
\r
- for(int i = 0;i<5;++i){\r
- sol += pow((x-gauss_n5[i]*y)*(x-gauss_n5[i]*y)+l*l,p)*gauss_w5[i]*y;\r
+ for (int i = 0; i < 5; ++i) {\r
+ sol += pow((x - gauss_n5[i] * y) * (x - gauss_n5[i] * y) + l * l, p)\r
+ * gauss_w5[i] * y;\r
}\r
\r
return sol;\r
double G_QY1Y2(double p, double y1, double y2, double x1, double x2, double l) {\r
\r
double sol = 0;\r
- for(int i = 0;i<5;++i){\r
- for(int j=0;j<5;++j){\r
- sol += pow((x1-y1*gauss_n5[i])*(x1-y1*gauss_n5[i])\r
- +(x2-y2*gauss_n5[j])*(x2-y2*gauss_n5[j])+l*l,p)*y1*gauss_w5[i]*y2*gauss_w5[j];\r
+ for (int i = 0; i < 5; ++i) {\r
+ for (int j = 0; j < 5; ++j) {\r
+ sol += pow(\r
+ (x1 - y1 * gauss_n5[i]) * (x1 - y1 * gauss_n5[i]) + (x2\r
+ - y2 * gauss_n5[j]) * (x2 - y2 * gauss_n5[j]) + l\r
+ * l, p) * y1 * gauss_w5[i] * y2 * gauss_w5[j];\r
}\r
}\r
return sol;\r
}\r
\r
-\r
double G_AY2X2(double y1, double y2, double x1, double x2, double d3) {\r
\r
// printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
\r
- double hlp = ((x1 - y1 ) * (x1 - y1 ) + (x2 - y2) * (x2 - y2) + d3 * d3);\r
+ double hlp = ((x1 - y1) * (x1 - y1) + (x2 - y2) * (x2 - y2) + d3 * d3);\r
\r
double sol = sqrt(hlp);\r
\r
if ((x2 - y2) != 0)\r
- sol += (x2 - y2) * log(sqrt(hlp)-(x2 - y2));\r
+ sol += (x2 - y2) * log(sqrt(hlp) - (x2 - y2));\r
\r
return sol;\r
}\r
sol /= 2 * p + 2;\r
} else {\r
sol = NAN;\r
- cout << "warning in G_AY1Y2: no case for p="<< p <<" defined. Possible: [-1.5,-0.5,0.5,...]\n";\r
+ cout << "warning in G_AY1Y2: no case for p=" << p\r
+ << " defined. Possible: [-1.5,-0.5,0.5,...]\n";\r
}\r
\r
return sol;\r
}\r
\r
-double Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2, double l) {\r
+double Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2,\r
+ double l) {\r
double sol = 0;\r
\r
//sol += y2*G_AY1Y2(y1,x2,)\r
return sol;\r
}\r
\r
-\r
-\r
double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,\r
double d3) {\r
\r
}\r
\r
/*double applyInt4(\r
- double(*f)(double, double, double, double, double, double, double),\r
- double a, double b, double c, double d, double r, double t, double u,\r
- double v, double d1, double d2, double d3) {\r
+ double(*f)(double, double, double, double, double, double, double),\r
+ double a, double b, double c, double d, double r, double t, double u,\r
+ double v, double d1, double d2, double d3) {\r
\r
- return f(b, d, t, v, d1, d2, d3) - f(b, d, t, u, d1, d2, d3) - f(b, d, r,\r
- v, d1, d2, d3) + f(b, d, r, u, d1, d2, d3) - f(b, c, t, v, d1, d2,\r
- d3) + f(b, c, t, u, d1, d2, d3) + f(b, c, r, v, d1, d2, d3) - f(b,\r
- c, r, u, d1, d2, d3) - f(a, d, t, v, d1, d2, d3) + f(a, d, t, u,\r
- d1, d2, d3) + f(a, d, r, v, d1, d2, d3) - f(a, d, r, u, d1, d2, d3)\r
- + f(a, c, t, v, d1, d2, d3) - f(a, c, t, u, d1, d2, d3) - f(a, c,\r
- r, v, d1, d2, d3) + f(a, c, r, u, d1, d2, d3);\r
+ return f(b, d, t, v, d1, d2, d3) - f(b, d, t, u, d1, d2, d3) - f(b, d, r,\r
+ v, d1, d2, d3) + f(b, d, r, u, d1, d2, d3) - f(b, c, t, v, d1, d2,\r
+ d3) + f(b, c, t, u, d1, d2, d3) + f(b, c, r, v, d1, d2, d3) - f(b,\r
+ c, r, u, d1, d2, d3) - f(a, d, t, v, d1, d2, d3) + f(a, d, t, u,\r
+ d1, d2, d3) + f(a, d, r, v, d1, d2, d3) - f(a, d, r, u, d1, d2, d3)\r
+ + f(a, c, t, v, d1, d2, d3) - f(a, c, t, u, d1, d2, d3) - f(a, c,\r
+ r, v, d1, d2, d3) + f(a, c, r, u, d1, d2, d3);\r
\r
-}*/\r
+ }*/\r
\r
double apply0Int4(\r
double(*f)(double, double, double, double, double, double, double),\r
double(*f)(double, double, double, double, double, double, double),\r
double b, double d, double t, double v, double d1, double d2, double d3) {\r
\r
- return f(b, d, t, v, d1, d2, d3)\r
- -f(b, 0, t, v, d1, d2, d3)\r
- -f(0, d, t, v, d1, d2, d3)\r
- +f(0, 0, t, v, d1, d2, d3);\r
+ return f(b, d, t, v, d1, d2, d3) - f(b, 0, t, v, d1, d2, d3) - f(0, d, t,\r
+ v, d1, d2, d3) + f(0, 0, t, v, d1, d2, d3);\r
\r
}\r
\r
// Berechnet das eigentliche Integral für parallele Elemente voll Analytisch\r
-double calcParIntA(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+double calcParIntA(double b, double d, double t, double v, double d1,\r
+ double d2, double d3) {\r
return apply0Int4(F_par, b, d, t, v, d1, d2, d3);\r
\r
}\r
\r
for (int i = 0; i < 5; ++i) {\r
for (int j = 0; j < 5; ++j) {\r
- sol += G_AY2X2( t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i],d , d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
- sol -= G_AY2X2( t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], 0, d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
- sol -= G_AY2X2( t * gauss_n5[j] + d1, d2, b * gauss_n5[i], d, d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
- sol += G_AY2X2( t * gauss_n5[j] + d1, d2, b * gauss_n5[i], 0, d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
+ sol\r
+ += G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i],\r
+ d, d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
+ sol\r
+ -= G_AY2X2(t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i],\r
+ 0, d3) * t * b * gauss_w5[i] * gauss_w5[j];\r
+ sol -= G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], d, d3)\r
+ * t * b * gauss_w5[i] * gauss_w5[j];\r
+ sol += G_AY2X2(t * gauss_n5[j] + d1, d2, b * gauss_n5[i], 0, d3)\r
+ * t * b * gauss_w5[i] * gauss_w5[j];\r
}\r
}\r
\r
\r
}\r
\r
-double calcParIntQY1(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+double calcParIntQY1(double b, double d, double t, double v, double d1,\r
+ double d2, double d3) {\r
//Fall 4\r
return apply0Int4(F_par, b, d, t, v, d1, d2, d3); ///ACHTUNG noch Falsch\r
\r
}\r
\r
-double calcOrtIntA(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+double calcOrtIntA(double b, double d, double t, double v, double d1,\r
+ double d2, double d3) {\r
return apply0Int4(F_ort, b, d, t, v, d1, d2, d3);\r
\r
}\r
\r
-double calcOrtIntQY1Y2(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+double calcOrtIntQY1Y2(double b, double d, double t, double v, double d1,\r
+ double d2, double d3) {\r
return 0;\r
\r
}\r
\r
+double calcParIntO0(double b, double d, double t, double v, double d1,\r
+ double d2, double d3, double eta) {\r
+ return calcParIntA(b,d,t,v,d1,d2,d3);\r
+}\r
+double calcOrtIntO0(double b, double d, double t, double v, double d1,\r
+ double d2, double d3, double eta) {\r
+ return calcOrtIntA(b,d,t,v,d1,d2,d3);\r
+}\r
\r
-double calcParIntO1(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+double calcParIntO1(double b, double d, double t, double v, double d1,\r
+ double d2, double d3, double eta) {\r
\r
//Elmente anordnen\r
- if(b*b+d*d>t*t+v*v){\r
+ if (b * b + d * d > t * t + v * v) {\r
double tmp = 0;\r
\r
- tmp = b; b = t; t = tmp;\r
- tmp = d; d = v; v = tmp;\r
+ tmp = b;\r
+ b = t;\r
+ t = tmp;\r
+ tmp = d;\r
+ d = v;\r
+ v = tmp;\r
\r
d1 = -d1;\r
d2 = -d2;\r
d3 = -d3;\r
}\r
\r
- return calcParIntA( b, d, t, v, d1, d2, d3);\r
+ return calcParIntA(b, d, t, v, d1, d2, d3);\r
}\r
\r
-double calcOrtIntO1(double b, double d, double t, double v, double d1, double d2, double d3) {\r
+double calcOrtIntO1(double b, double d, double t, double v, double d1,\r
+ double d2, double d3, double eta) {\r
\r
//Elmente anordnen\r
- if(b*b+d*d>t*t+v*v){\r
+ if (b * b + d * d > t * t + v * v) {\r
double tmp = 0;\r
\r
- tmp = b; b = t; t = tmp;\r
- tmp = d; d = v; v = tmp;\r
+ tmp = b;\r
+ b = t;\r
+ t = tmp;\r
+ tmp = d;\r
+ d = v;\r
+ v = tmp;\r
\r
d1 = -d1;\r
d2 = -d2;\r
d3 = -d3;\r
}\r
\r
- return calcOrtIntA( b, d, t, v, d1, d2, d3);\r
+ return calcOrtIntA(b, d, t, v, d1, d2, d3);\r
}\r
\r
-\r
-double calcParIntO2(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {\r
+double calcParIntO2(double b, double d, double t, double v, double d1,\r
+ double d2, double d3, double eta) {\r
\r
//Elmente anordnen\r
- if(b*b+d*d>t*t+v*v){\r
+ if (b * b + d * d > t * t + v * v) {\r
double tmp = 0;\r
\r
- tmp = b; b = t; t = tmp;\r
- tmp = d; d = v; v = tmp;\r
+ tmp = b;\r
+ b = t;\r
+ t = tmp;\r
+ tmp = d;\r
+ d = v;\r
+ v = tmp;\r
\r
d1 = -d1;\r
d2 = -d2;\r
d3 = -d3;\r
}\r
\r
- if((b*b+d*d) *eta <= d1*d1+d2*d2+d3*d3){\r
- return calcParIntQX1X2( b, d, t, v, d1, d2, d3);\r
- }else{\r
- return calcParIntA( b, d, t, v, d1, d2, d3);\r
+ if ((b * b + d * d) * eta <= d1 * d1 + d2 * d2 + d3 * d3) {\r
+ return calcParIntQX1X2(b, d, t, v, d1, d2, d3);\r
+ } else {\r
+ return calcParIntA(b, d, t, v, d1, d2, d3);\r
}\r
\r
-\r
}\r
\r
-double calcOrtIntO2(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {\r
+double calcOrtIntO2(double b, double d, double t, double v, double d1,\r
+ double d2, double d3, double eta) {\r
\r
//Elmente anordnen\r
- if(b*b+d*d>t*t+v*v){\r
+ if (b * b + d * d > t * t + v * v) {\r
double tmp = 0;\r
\r
- tmp = b; b = t; t = tmp;\r
- tmp = d; d = v; v = tmp;\r
+ tmp = b;\r
+ b = t;\r
+ t = tmp;\r
+ tmp = d;\r
+ d = v;\r
+ v = tmp;\r
\r
d1 = -d1;\r
d2 = -d2;\r
d3 = -d3;\r
}\r
\r
- if((b*b+d*d) *eta <= d1*d1+d2*d2+d3*d3){\r
+ if ((b * b + d * d) * eta <= d1 * d1 + d2 * d2 + d3 * d3) {\r
return 0;\r
- }else{\r
- return calcOrtIntA( b, d, t, v, d1, d2, d3);\r
+ } else {\r
+ return calcOrtIntA(b, d, t, v, d1, d2, d3);\r
}\r
}\r
\r
-\r
-double calcParIntO3(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {\r
+double calcParIntO3(double b, double d, double t, double v, double d1,\r
+ double d2, double d3, double eta) {\r
\r
//Achsen vertauschen\r
- if(b*b+t*t>d*d+v*v){\r
+ if (b * b + t * t > d * d + v * v) {\r
double tmp = 0;\r
\r
- tmp = b; b = d; d = tmp;\r
- tmp = t; t = v; v = tmp;\r
+ tmp = b;\r
+ b = d;\r
+ d = tmp;\r
+ tmp = t;\r
+ t = v;\r
+ v = tmp;\r
\r
- tmp = d1; d1 = d2; d2 = tmp;\r
+ tmp = d1;\r
+ d1 = d2;\r
+ d2 = tmp;\r
}\r
\r
- if(max(b,t) *eta <= d1){\r
+ if (max(b, t) * eta <= d1) {\r
return calcParIntQY1X1(b, d, t, v, d1, d2, d3);\r
- }else{\r
- return calcParIntA( b, d, t, v, d1, d2, d3);\r
+ } else {\r
+ return calcParIntA(b, d, t, v, d1, d2, d3);\r
}\r
\r
-\r
}\r
\r
-double calcOrtIntO3(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {\r
+double calcOrtIntO3(double b, double d, double t, double v, double d1,\r
+ double d2, double d3, double eta) {\r
\r
//Elmente anordnen\r
- if(b*b+d*d>t*t+v*v){\r
+ if (b * b + d * d > t * t + v * v) {\r
double tmp = 0;\r
\r
- tmp = b; b = t; t = tmp;\r
- tmp = d; d = v; v = tmp;\r
+ tmp = b;\r
+ b = t;\r
+ t = tmp;\r
+ tmp = d;\r
+ d = v;\r
+ v = tmp;\r
\r
d1 = -d1;\r
d2 = -d2;\r
d3 = -d3;\r
}\r
\r
- if(max(b,t) *eta > d1){\r
+ if (max(b, t) * eta > d1) {\r
return 0;\r
- }else{\r
- return calcOrtIntA( b, d, t, v, d1, d2, d3);\r
+ } else {\r
+ return calcOrtIntA(b, d, t, v, d1, d2, d3);\r
}\r
}\r
#ifndef HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
#define HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
-//int sign(double);
-//double arsinh(double);
-
-
// sol = g0(p,y,x,l);
double g_AY(double, double, double, double);
// sol = g0(p,y,x,l);
double g_QY(double, double, double, double);
-
// sol = G00(p,y1,y2,x1,x2,l);
double G_AY1Y2(double, double, double, double, double, double);
// sol = G00(p,y1,y2,x1,x2,l);
// sol = G00(p,y1,y2,x1,x2,l);
double G_QY1Y2(double, double, double, double, double, double);
-
// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
double F_par(double, double, double, double, double, double, double);
// sol = F_ort(x1,x2,y2,y3,d1,d2,d3);
double F_ort(double, double, double, double, double, double, double);
-
// sol = quad0Int4((F_par/F_ort),b,d,t,v,d1,d2,d3);
double apply0Int4(
double(*)(double, double, double, double, double, double, double),
double calcParIntA(double, double, double, double, double, double, double);
double calcParIntQX1X2(double, double, double, double, double, double, double);
double calcParIntQY1X1(double, double, double, double, double, double, double);
-
double calcParIntQX1(double, double, double, double, double, double, double);
-
// sol = calcOrtInt.(b,d,t,v,d1,d2,d3);
double calcOrtIntA(double, double, double, double, double, double, double);
double calcOrtIntQX1X2(double, double, double, double, double, double, double);
+//Voll Analytische Berechnung
+double calcParIntO0(double, double, double, double, double, double, double,
+ double);
+double calcOrtIntO0(double, double, double, double, double, double, double,
+ double);
+
// Elemente Vertauschen
// sol = calc...Int01(b,d,t,v,d1,d2,d3);
-double calcParIntO1(double, double, double, double, double, double, double);
-double calcOrtIntO1(double, double, double, double, double, double, double);
+double calcParIntO1(double, double, double, double, double, double, double,
+ double);
+double calcOrtIntO1(double, double, double, double, double, double, double,
+ double);
// Quadratur ueber kleineres Element
// sol = calc...Int02(b,d,t,v,d1,d2,d3,eta);
-double calcParIntO2(double, double, double, double, double, double, double, double);
-double calcOrtIntO2(double, double, double, double, double, double, double, double);
+double calcParIntO2(double, double, double, double, double, double, double,
+ double);
+double calcOrtIntO2(double, double, double, double, double, double, double,
+ double);
// Quadratur kuerzere Achse/Seite
// sol = calc...Int03(b,d,t,v,d1,d2,d3,eta);
-double calcParIntO3(double, double, double, double, double, double, double, double);
-double calcOrtIntO3(double, double, double, double, double, double, double, double);
+double calcParIntO3(double, double, double, double, double, double, double,
+ double);
+double calcOrtIntO3(double, double, double, double, double, double, double,
+ double);
#endif
+++ /dev/null
-function [ret] = test_Fort(a,b,c,d,r,t,u,v,l1,l2,l3)
-
-
- quad_sol = surfDoubleQuad(@(x1,x2,y2,y3)1/sqrt((x1-l1).^2+(x2-y2-l2).^2+(y3+l3).^2)...
- ,a,b,c,d,r,t,u,v,4);
-
- 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);
-
- sol = intF(@(x1,x2,y1,y2)mex_Fort(x1,x2,y1,y2,l1,l2,l3),a,b,c,d,r,t,u,v);
-
-
-
- ret = [quad_sol sol];
-end
-
+++ /dev/null
-function [ret] = test_Fpar(a,b,c,d,r,t,u,v,l1,l2,l3)
-%
-% [ret] = test_Fpar(a,b,c,d,r,t,u,v,l1,l2,l3)
-
- quad_sol = surfDoubleQuad(@(x1,x2,y1,y2)1/sqrt((x1-y1-l1).^2+(x2-y2-l2).^2+l3.^2)...
- ,a,b,c,d,r,t,u,v,5);
-
- 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);
-
-
-
-
- [ mex_Fpar(b,d,t,v,l1,l2,l3) mex_FparY1Y2_2(b,d,t,v,l1,l2,l3);...
- mex_Fpar(b,d,t,u,l1,l2,l3) mex_FparY1Y2_2(b,d,t,u,l1,l2,l3);...
- mex_Fpar(b,d,r,v,l1,l2,l3) mex_FparY1Y2_2(b,d,r,v,l1,l2,l3)];
-
- [i j] = mex_G00(-0.5, b, d, t + l1, v + l2, l3);
-
- sol = intF(@(x1,x2,y1,y2)mex_Fpar(x1,x2,y1,y2,l1,l2,l3),a,b,c,d,r,t,u,v);
-
- sol2 = intF(@(x1,x2,y1,y2)mex_FparY1Y2_2(x1,x2,y1,y2,l1,l2,l3),a,b,c,d,r,t,u,v);
-
- int2 = @(x1,y1,c,d,u,v)...
- mex_FparQY2X2(x1,d,y1,v,l1,l2,l3)...
- -mex_FparQY2X2(x1,c,y1,v,l1,l2,l3)...
- -mex_FparQY2X2(x1,d,y1,u,l1,l2,l3)...
- +mex_FparQY2X2(x1,c,y1,u,l1,l2,l3);
-
- sol3 = surfQuad(@(x1,y1)(int2(x1,y1,c,d,u,v)),a,b,r,t,5);
-
- sol4 =...
- surfQuad(@(x1,y1)mex_FparQY2X2(x1,d,y1,v,l1,l2,l3),a,b,r,t,5)...
- -surfQuad(@(x1,y1)mex_FparQY2X2(x1,c,y1,v,l1,l2,l3),a,b,r,t,5)...
- -surfQuad(@(x1,y1)mex_FparQY2X2(x1,d,y1,u,l1,l2,l3),a,b,r,t,5)...
- +surfQuad(@(x1,y1)mex_FparQY2X2(x1,c,y1,u,l1,l2,l3),a,b,r,t,5);
-
-
- ret = [quad_sol sol sol2 sol3 sol4];
-end
-
+++ /dev/null
-function [ret] = test_G00(p, x1, x2,l,a,b,c,d)
-%
-% [ret] = test_G00(p, x1, x2,l,a,b,c,d)
-
- quad_sol = surfQuad(@(y1,y2)((x1-y1).^2+(x2-y2).^2+l.^2).^p,a,b,c,d,4);
-
-
- [my(1) mq(1)] = mex_G00(p,a,c,x1,x2,l)
- [my(2) mq(2)] = mex_G00(p,b,d,x1,x2,l)
- [my(3) mq(3)] = mex_G00(p,b,c,x1,x2,l)
- [my(4) mq(4)] = mex_G00(p,a,d,x1,x2,l)
-
-
- my_sol = my(1)+my(2)-my(3)-my(4);
- mq_sol = mq(1)+mq(2)-mq(3)-mq(4);
-
- ret = [quad_sol my_sol mq_sol];
-end
-
+++ /dev/null
-function done=test_functions(eps)
-
- done = 0;
-%% g0 testen
-p_g0 = [0.5 0 -0.5 -1 -1.5];
-for l= 1:5
- for p = p_g0
- s = test_g0(p,1.4,l,1,4);
- if(abs(s(1)-s(2))>eps)
- disp g0
- l
- p
- s
- return;
- end
- end
-end
- % x muss ausserhalb des Intervalls liegen
- for p = p_g0
- s = test_g0(p,1.4,0,2,4);
- if(abs(s(1)-s(2))>eps)
- disp g0
- p
- s
- return;
- end
- end
-
-%% G00 testen
-p_G00 = [-0.5 -1.5];
-
- for p = p_G00
- s = test_G00(p,1.4,0.7,1,2,4,1,3);
- if(abs(s(1)-s(2))>eps)
- disp G00
- p
- s
- return;
- end
- end
-
- for p = p_G00
- s = test_G00(p,1.1,0.3,0,2,4.5,1,2.5);
- if(abs(s(1)-s(2))>eps)
- disp G00
- p
- s
- return;
- end
- end
-
-%% F_par testen
-
-
- s = test_Fpar(1,3,2,3,4,5,1,3,1,3,1);
- if(abs(s(1)-s(2))>eps)
- disp Fpar
- p
- s
- return;
- end
-
-
-
- done = 1;
-end
\ No newline at end of file
+++ /dev/null
-function [ret] = test_g0(p, x,l,a,b)
-
-
-% quad_sol = quad(@(y)((x-y).^2+l.^2).^p,a,b);
-[no we] = gauss(5,a,b);
-f = @(y)((x-y).^2+l.^2).^p;
-quad_sol = we*f(no);
-
-% cumsum(f(no))'
-
-% X = 0:0.2:5;
-% Y = ((x-X).^2+l.^2).^p;
-%
-% plot(X,Y)
-
-[my(1) mq(1)] = mex_g0(p,a,x,l);
-[my(2) mq(2)]= mex_g0(p,b,x,l);
-
-my_sol = my(2) - my(1);
-mq_sol = mq(2) - mq(1);
-
-
-
-ret = [quad_sol my_sol mq_sol];
-end
-
+++ /dev/null
-\r
-mex mex_build_As1.cpp slpRectangle.cpp\r
-\r
-dataIso =[];\r
-dataAniso =[];\r
-maxSize = 10^3;\r
-Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren\r
-% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren\r
-% Netz = 'exmpl_2DQuad';\r
-\r
-%% Isotrope Verfeinerung\r
-disp 'Isotrop'\r
-load(Netz)\r
-\r
-ref = 0;\r
-while size(elements,1)<maxSize\r
- tic\r
- ref = ref+1;\r
-\r
- [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
- A_fine = mex_build_As1(coordinates_fine,elements_fine);\r
- b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
- x_fine = A_fine\b;\r
- xe_fine = x_fine'*A_fine*x_fine;\r
-\r
- ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
- marked = 2*ones(1,size(elements,1));\r
-\r
- dataIso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
- [coordinates, elements] = refineQuad(coordinates,elements,marked);\r
- toc\r
-end\r
-dataIso\r
-\r
-\r
-%% Anisotrope Verfeinerung\r
-disp 'Anisotrop'\r
-load(Netz)\r
-\r
-ref = 0;\r
-while size(elements,1)<maxSize\r
- tic\r
- ref = ref+1;\r
-\r
- %Netz Isotrop verfeinern (4 Teile)\r
- [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
- \r
- %Matrix mit feinem Netz berechnen\r
- A_fine = mex_build_As1(coordinates_fine,elements_fine);\r
-% A2 = mex_build_A(coordinates_fine,elements_fine);\r
- \r
-% figure(1)\r
-% spy(A_fine-A_fine')\r
-% \r
-% figure (2)\r
-% spy(A2-A2')\r
-% eh = A_fine - A2;\r
-% sum(sum(abs(A2-A2')))\r
-% sum(sum(abs(A_fine-A_fine')))\r
- \r
-% spy(eh)\r
- \r
-% sum(sum(abs(eh)))\r
-% if(eh.^2)\r
-% disp 'hä'\r
-% end\r
- \r
- %rechte Seite Aufbauen\r
- b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
- \r
- %Lösungsvektor bestimmen\r
- x_fine = A_fine\b;\r
- \r
- %Energienorm^2 bestimmen\r
- xe_fine = x_fine'*A_fine*x_fine;\r
-\r
- % Welche Elemente sollen verfeinert werden\r
- ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
- \r
- % Wie soll verfeinert werden\r
- marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
-\r
- %Daten zur Auswertung zwischenspeichern (testweise?)\r
- dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
- \r
- %Sicherheit falls keine Elemente markiert sind\r
- if(sum(marked~=1)==0)\r
- disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
- break;\r
- end\r
- \r
- %Verfeinern\r
- [coordinates, elements] = refineQuad(coordinates,elements,marked);\r
- toc\r
-end\r
-dataAniso\r
-\r
-clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine\r
-\r
-\r
-%% Alles Zeichnen\r
-figure (1)\r
-loglog(dataIso(:,1),dataIso(:,2),dataAniso(:,1),dataAniso(:,2));\r
-\r
-legend('Isotrop','Anisotrop')\r
-xlabel('Elemente')\r
-ylabel('Fehler')\r
-\r
-figure (2)\r
-loglog(dataIso(:,1),dataIso(:,3),dataAniso(:,1),dataAniso(:,3));\r
-\r
-legend('Isotrop','Anisotrop')\r
-xlabel('Elemente')\r
-ylabel('EnergieNorm')\r
+++ /dev/null
-\r
-mex mex_build_As2.cpp slpRectangle.cpp\r
-\r
-dataIso =[];\r
-dataAniso =[];\r
-maxSize = 500;\r
-fileName = './test_save';\r
-Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren\r
-% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren\r
-% Netz = 'exmpl_2DQuad';\r
-\r
-\r
-%% Anisotrope Verfeinerung\r
-\r
-if(~exist('d1'))\r
- \r
- disp 'Normal'\r
- load(Netz)\r
-\r
- ref = 0;\r
-\r
-while size(elements,1)<maxSize\r
- tic\r
- ref = ref+1; \r
- \r
-\r
-\r
- %Netz Isotrop verfeinern (4 Teile)\r
- [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
- \r
- %Matrix mit feinem Netz berechnen\r
- A_fine = mex_build_A(coordinates_fine,elements_fine);\r
- \r
- %rechte Seite Aufbauen\r
- b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
- \r
- %Lösungsvektor bestimmen\r
- x_fine = A_fine\b;\r
- \r
- %Energienorm^2 bestimmen\r
- xe_fine = x_fine'*A_fine*x_fine;\r
-\r
- % Welche Elemente sollen verfeinert werden\r
- ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
- \r
- % Wie soll verfeinert werden\r
- marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
-\r
- %Daten zur Auswertung zwischenspeichern (testweise?)\r
- dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
- \r
- %Sicherheit falls keine Elemente markiert sind\r
- if(sum(marked~=1)==0)\r
- disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
- break;\r
- end\r
- \r
- %Verfeinern\r
- %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
- file = ['save' int2str(ref)];\r
- load(file);\r
- toc\r
-end\r
-d1 = dataAniso\r
-end\r
-\r
-%% Anisotrope Verfeinerung\r
-disp 'Optimiert'\r
-load(Netz)\r
-\r
-ref = 0;\r
-while size(elements,1)<maxSize\r
- tic\r
- ref = ref+1; \r
- \r
-\r
-\r
- %Netz Isotrop verfeinern (4 Teile)\r
- [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
- \r
- %Matrix mit feinem Netz berechnen\r
- A_fine = mex_build_As2(coordinates_fine,elements_fine);\r
- \r
- %rechte Seite Aufbauen\r
- b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
- \r
- %Lösungsvektor bestimmen\r
- x_fine = A_fine\b;\r
- \r
- %Energienorm^2 bestimmen\r
- xe_fine = x_fine'*A_fine*x_fine;\r
-\r
- % Welche Elemente sollen verfeinert werden\r
- ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
- \r
- % Wie soll verfeinert werden\r
- marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
-\r
- %Daten zur Auswertung zwischenspeichern (testweise?)\r
- dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
- \r
- %Sicherheit falls keine Elemente markiert sind\r
-% if(sum(marked~=1)==0)\r
-% disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
-% break;\r
-% end\r
- \r
- %Verfeinern\r
- %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
- file = ['save' int2str(ref)];\r
- load(file);\r
- toc\r
-end\r
-d2 = dataAniso\r
-\r
-clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine\r
-\r
-\r
-%% Alles Zeichnen\r
-figure (1)\r
-loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));\r
-\r
-legend('Normal','optimiert')\r
-xlabel('Elemente')\r
-ylabel('Fehler')\r
-\r
-figure (2)\r
-loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));\r
-\r
-legend('Normal','Optimiert')\r
-xlabel('Elemente')\r
-ylabel('EnergieNorm')\r
+++ /dev/null
-\r
-mex mex_build_As3.cpp slpRectangle.cpp\r
-\r
-dataIso =[];\r
-dataAniso =[];\r
-maxSize = 500;\r
-fileName = './test_save';\r
-Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren\r
-% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren\r
-% Netz = 'exmpl_2DQuad';\r
-\r
-\r
-%% Anisotrope Verfeinerung\r
-\r
-if(~exist('d1'))\r
- \r
- disp 'Normal'\r
- load(Netz)\r
-\r
- ref = 0;\r
-\r
-while size(elements,1)<maxSize\r
- tic\r
- ref = ref+1; \r
- \r
-\r
-\r
- %Netz Isotrop verfeinern (4 Teile)\r
- [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
- \r
- %Matrix mit feinem Netz berechnen\r
- A_fine = mex_build_A(coordinates_fine,elements_fine);\r
- \r
- %rechte Seite Aufbauen\r
- b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
- \r
- %Lösungsvektor bestimmen\r
- x_fine = A_fine\b;\r
- \r
- %Energienorm^2 bestimmen\r
- xe_fine = x_fine'*A_fine*x_fine;\r
-\r
- % Welche Elemente sollen verfeinert werden\r
- ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
- \r
- % Wie soll verfeinert werden\r
- marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
-\r
- %Daten zur Auswertung zwischenspeichern (testweise?)\r
- dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
- \r
- %Sicherheit falls keine Elemente markiert sind\r
- if(sum(marked~=1)==0)\r
- disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
- break;\r
- end\r
- \r
- %Verfeinern\r
- %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
- file = ['save' int2str(ref)];\r
- load(file);\r
- toc\r
-end\r
-d1 = dataAniso\r
-end\r
-\r
-%% Anisotrope Verfeinerung\r
-disp 'Optimiert'\r
-load(Netz)\r
-\r
-ref = 0;\r
-while size(elements,1)<maxSize\r
- tic\r
- ref = ref+1; \r
- \r
-\r
-\r
- %Netz Isotrop verfeinern (4 Teile)\r
- [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
- \r
- %Matrix mit feinem Netz berechnen\r
- A_fine = mex_build_As3(coordinates_fine,elements_fine);\r
- \r
- %rechte Seite Aufbauen\r
- b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
- \r
- %Lösungsvektor bestimmen\r
- x_fine = A_fine\b;\r
- \r
- %Energienorm^2 bestimmen\r
- xe_fine = x_fine'*A_fine*x_fine;\r
-\r
- % Welche Elemente sollen verfeinert werden\r
- ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
- \r
- % Wie soll verfeinert werden\r
- marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
-\r
- %Daten zur Auswertung zwischenspeichern (testweise?)\r
- dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
- \r
- %Sicherheit falls keine Elemente markiert sind\r
-% if(sum(marked~=1)==0)\r
-% disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
-% break;\r
-% end\r
- \r
- %Verfeinern\r
- %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
- file = ['save' int2str(ref)];\r
- load(file);\r
- toc\r
-end\r
-d2 = dataAniso\r
-\r
-clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine\r
-\r
-\r
-%% Alles Zeichnen\r
-figure (1)\r
-loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));\r
-\r
-legend('Normal','optimiert')\r
-xlabel('Elemente')\r
-ylabel('Fehler')\r
-\r
-figure (2)\r
-loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));\r
-\r
-legend('Normal','Optimiert')\r
-xlabel('Elemente')\r
-ylabel('EnergieNorm')\r
--- /dev/null
+\r
+mex mex_build_AU.cpp slpRectangle.cpp\r
+\r
+dataIso =[];\r
+dataAniso =[];\r
+maxSize = 500;\r
+fileName = './test_save';\r
+Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren\r
+% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 16.2265 konvergieren\r
+% Netz = 'exmpl_2DQuad';\r
+\r
+\r
+%% Anisotrope Verfeinerung\r
+\r
+if(~exist('d1'))\r
+ \r
+ disp 'Normal'\r
+ load(Netz)\r
+\r
+ ref = 0;\r
+\r
+while size(elements,1)<maxSize\r
+ tic\r
+ ref = ref+1; \r
+ \r
+ %Netz Isotrop verfeinern (4 Teile)\r
+ [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
+ \r
+ %Matrix mit feinem Netz berechnen\r
+ A_fine = mex_build_A(coordinates_fine,elements_fine);\r
+ \r
+ %rechte Seite Aufbauen\r
+ b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
+ \r
+ %Lösungsvektor bestimmen\r
+ x_fine = A_fine\b;\r
+ \r
+ %Energienorm^2 bestimmen\r
+ xe_fine = x_fine'*A_fine*x_fine;\r
+\r
+ % Welche Elemente sollen verfeinert werden\r
+ ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
+ \r
+ % Wie soll verfeinert werden\r
+ marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
+\r
+ %Daten zur Auswertung zwischenspeichern (testweise?)\r
+ dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
+ \r
+ %Sicherheit falls keine Elemente markiert sind\r
+ if(sum(marked~=1)==0)\r
+ disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
+ break;\r
+ end\r
+ \r
+ %Verfeinern\r
+ %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
+ file = ['save' int2str(ref)];\r
+ load(file);\r
+ toc\r
+end\r
+d1 = dataAniso\r
+end\r
+\r
+%% Anisotrope Verfeinerung\r
+disp 'Optimiert'\r
+load(Netz)\r
+\r
+ref = 0;\r
+while size(elements,1)<maxSize\r
+ tic\r
+ ref = ref+1; \r
+ \r
+\r
+\r
+ %Netz Isotrop verfeinern (4 Teile)\r
+ [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
+ \r
+ %Matrix mit feinem Netz berechnen\r
+ A_fine = mex_build_AU(coordinates_fine,elements_fine);\r
+ \r
+ %rechte Seite Aufbauen\r
+ b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));\r
+ \r
+ %Lösungsvektor bestimmen\r
+ x_fine = A_fine\b;\r
+ \r
+ %Energienorm^2 bestimmen\r
+ xe_fine = x_fine'*A_fine*x_fine;\r
+\r
+ % Welche Elemente sollen verfeinert werden\r
+ ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);\r
+ \r
+ % Wie soll verfeinert werden\r
+ marked = mark(x_fine(f2s)',ind,0.5,0.5);\r
+\r
+ %Daten zur Auswertung zwischenspeichern (testweise?)\r
+ dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];\r
+ \r
+ %Sicherheit falls keine Elemente markiert sind\r
+% if(sum(marked~=1)==0)\r
+% disp 'Abgebrochen weil nicht mehr verfeinert wird!'\r
+% break;\r
+% end\r
+ \r
+ %Verfeinern\r
+ %[coordinates, elements] = refineQuad(coordinates,elements,marked);\r
+ file = ['save' int2str(ref)];\r
+ load(file);\r
+ toc\r
+end\r
+d2 = dataAniso\r
+\r
+clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine\r
+\r
+\r
+%% Alles Zeichnen\r
+figure (1)\r
+loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));\r
+\r
+legend('Normal','optimiert')\r
+xlabel('Elemente')\r
+ylabel('Fehler')\r
+\r
+figure (2)\r
+loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));\r
+\r
+legend('Normal','Optimiert')\r
+xlabel('Elemente')\r
+ylabel('EnergieNorm')\r
+++ /dev/null
-
-N = size(elements,1);
-
-A1 = zeros(N);
-
-x = zeros(N,2);
-xe = zeros(1,2);
-t = zeros(1,2);
-
-% untere schranke s t obere schranke k l
-intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ...
- f(k1,k2,l1,l2)-f(k1,k2,l1,t2)-f(k1,k2,t1,l2)+f(k1,k2,t1,t2)...
- -f(k1,s2,l1,l2)+f(k1,s2,l1,t2)+f(k1,s2,t1,l2)-f(k1,s2,t1,t2)...
- -f(s1,k2,l1,l2)+f(s1,k2,l1,t2)+f(s1,k2,t1,l2)-f(s1,l2,t1,t2)...
- +f(s1,s2,l1,l2)-f(s1,0,l1,t2)-f(s1,s2,t1,l2)+f(s1,s2,t1,t2);
-tic
-A2 = build_A2(coordinates,elements);
-t(1) = toc;
-disp ' '
-tic
-A2 = mex_build_A(coordinates,elements);
-t(2) = toc;
-b = sqrt(sum(quadNorm(coordinates,elements,'w').^2,2));
-
-x(:,1) = A1\b;
-xe(1) = x(:,1)'*A1*x(:,1);
-
-x(:,2) = A2\b;
-xe(2) = x(:,2)'*A2*x(:,2);
-
-x;
-xe
-t
\ No newline at end of file
+++ /dev/null
-\r
-rt_times = 5;\r
-\r
-time = zeros(rt_times,3);\r
-\r
- test_solveS\r
- \r
- time(1,:) = [size(x,1) t];\r
-\r
-for l = 2:rt_times\r
- test_refine\r
- test_solveS\r
- \r
- time(l,:) = [size(x,1) t];\r
-end\r
-\r
-\r
-loglog(time(:,1),time(:,2),time(:,1),time(:,3));\r
-\r
-legend('Matlab','Mex')\r
-xlabel('Elemente')\r
-ylabel('Zeit')\r
-\r
projects / bacc.git / commitdiff