int sign(double);\r
//double arsinh(double);\r
\r
-/*\r
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
- //sicherheitsabfragen zu Datengroessen\r
- if (nrhs != 2)\r
- mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx3))");\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 != 3)\r
- mexErrMsgTxt("expected elements (Mx3)");\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 * 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
-\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
-\r
- double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double tmp;\r
-\r
- int rx, ry;\r
-\r
- //Ausrechnen\r
- for (int j = 0; j < em; ++j) {\r
- x1[0] = C[(int) E[j] - 1];\r
- x1[1] = C[cm + (int) E[j] - 1];\r
- x1[2] = C[2 * cm + (int) E[j] - 1];\r
-\r
- x2[0] = C[(int) E[em + j] - 1];\r
- x2[1] = C[cm + (int) E[em + j] - 1];\r
- x2[2] = C[2 * cm + (int) E[em + j] - 1];\r
-\r
- x3[0] = C[(int) E[2 * em + j] - 1];\r
- x3[1] = C[cm + (int) E[2 * em + j] - 1];\r
- x3[2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
-\r
- xn[0] = (x2[1] - x1[1]) * (x3[2] - x1[2]) - (x2[2] - x1[2]) * (x3[1]\r
- - x1[1]);\r
- xn[1] = (x2[2] - x1[2]) * (x3[0] - x1[0]) - (x2[0] - x1[0]) * (x3[2]\r
- - x1[2]);\r
- xn[2] = (x2[0] - x1[0]) * (x3[1] - x1[1]) - (x2[1] - x1[1]) * (x3[0]\r
- - x1[0]);\r
-\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
- // printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);\r
-\r
- for (int k = 0; k < em; ++k) {\r
- y1[0] = C[(int) E[k] - 1];\r
- y1[1] = C[cm + (int) E[k] - 1];\r
- y1[2] = C[2 * cm + (int) E[k] - 1];\r
-\r
- y2[0] = C[(int) E[em + k] - 1];\r
- y2[1] = C[cm + (int) E[em + k] - 1];\r
- y2[2] = C[2 * cm + (int) E[em + k] - 1];\r
-\r
- y3[0] = C[(int) E[2 * em + k] - 1];\r
- y3[1] = C[cm + (int) E[2 * em + k] - 1];\r
- y3[2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
-\r
- yn[0] = (y2[1] - y1[1]) * (y3[2] - y1[2]) - (y2[2] - y1[2])\r
- * (y3[1] - y1[1]);\r
- yn[1] = (y2[2] - y1[2]) * (y3[0] - y1[0]) - (y2[0] - y1[0])\r
- * (y3[2] - y1[2]);\r
- yn[2] = (y2[0] - y1[0]) * (y3[1] - y1[1]) - (y2[1] - y1[1])\r
- * (y3[0] - y1[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
- d[0] = y1[0] - x1[0];\r
- d[1] = y1[1] - x1[1];\r
- d[2] = y1[2] - x1[2];\r
-\r
- if (rx = ry) {\r
- if (rx = 2) {\r
- //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x1[0],x1[1],x2[0],x3[1],y1[0],y1[1],y2[0],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 = quadInt(FLO_plane, x1[0], x1[1], x2[0], x3[1], y1[0],\r
- y1[1], y2[0], y3[1], d[0], d[1], d[2]);\r
- // printf("%d%d|%.2f\n",j,k,tmp);\r
- A[(k * em) + j] = 1. / (4 * my_PI) * tmp;\r
-\r
- } else\r
- A[(k * em) + j] = 0;\r
- } else {\r
- A[(k * em) + j] = 0;\r
- }\r
-\r
- // Vorbereiten der DATEN\r
- // ej = coordinates(elements(j,:)',:);\r
- // ek = coordinates(elements(k,:)',:);\r
- // d = ek(1,:) - ej(1,:);\r
-\r
- // AUSRECHNEN\r
- // A[(k*em)+j] = 1./(4*my_PI);// *\r
- // quadInt(F_par,\r
- // ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));\r
-\r
- }\r
-\r
- }\r
- //Rueckgabe (eventuell zurueck schreiben)\r
- // mxFree(x1);\r
- //mxFree(x2);\r
- //mxFree(x3);\r
- //mxFree(y1);\r
- //mxFree(y2);\r
- //mxFree(y3);\r
- //mxFree(xn);\r
- //mxFree(yn);\r
- //mxFree(d);\r
-\r
- return;\r
-}\r
-*/\r
int inline sign(double x) {\r
return x > 0 ? 1 : (x < 0 ? -1 : 0);\r
}\r
/*\r
-double inline arsinh(double x) {\r
- return log(x + sqrt(x * x + 1));\r
-}\r
-*/\r
+ double inline arsinh(double x) {\r
+ return log(x + sqrt(x * x + 1));\r
+ }\r
+ */\r
//y-x muss != 0 sein\r
double g0(double p, double y, double x, double l) {\r
//printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);\r
else if (p == -0.5)\r
sol = asinh((y - x) / fabs(l));\r
else if (p == -1)\r
- sol = atan((y - x) / fabs(l))/l;\r
+ sol = atan((y - x) / fabs(l)) / l;\r
else if (p == -1.5)\r
sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l));\r
else\r
\r
double G00(double p, double y1, double y2, double x1, double x2, double l) {\r
// printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);\r
-\r
+ double pt = p + 1.5;\r
double sol = 0;\r
- if (p == -1.5) {\r
+ if (pt == 0) {\r
if (l == 0) {\r
sol = -sqrt((y1 - x1) * (y1 - x1) + (y2 - x2) * (y2 - x2)) / ((y1\r
- x1) * (y2 - x2));\r
+ l * l)) + 1);\r
}\r
\r
- } else if (p == -0.5) {\r
+ } else if ((pt > 0) && ((int) pt == pt)) {\r
if (l != 0)\r
sol = 2 * p * l * l * G00(p - 1, y1, y2, x1, x2, l);\r
if ((y1 - x1) != 0)\r
sol /= 2 * p + 2;\r
} else {\r
sol = NAN;\r
- printf("warning in G00: no case for p defined\n");\r
+ printf("warning in G00: no case for p=%.2f defined\n", p);\r
//mexErrMsgTxt("no case for p defined\n");\r
}\r
\r
double sol = (x1 - y1 - d1) * (x2 - y2 - d2);\r
\r
if (sol != 0)\r
- sol *= G00(-0.5,x1, x2, y1 + d1, y2 + d2, d3);\r
+ sol *= G00(-0.5, x1, x2, y1 + d1, y2 + d2, d3);\r
\r
if ((x1 - y1 - d1) != 0)\r
sol -= (x1 - y1 - d1) * g0(0.5, x1, y1 + d1,\r
return sol;\r
}\r
\r
+double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,\r
+ double d3) {\r
+\r
+ // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y2,y3,d1,d2,d3);\r
+ double sol = -G00(0.5, y3, x1, -d3, d1, x2 - y2 - d2);\r
+\r
+ if ((x1 - d1) * (x2 - y2 - d2) != 0)\r
+ sol -= (x1 - d1) * (x2 - y2 - d2) * G00(-0.5, x2, y3, y2 + d2, -d3,\r
+ x1 - d1);\r
+\r
+ if ((x1 - d1) != 0)\r
+ sol += (x1 - d1) * g0(0.5, y3, -d3,\r
+ sqrt((x1 - d1) * (x1 - d1) + (x2 - y2 - d2) * (x2 - y2 - d2)));\r
+\r
+ if ((y3 + d3) * (x2 - y2 - d2) != 0)\r
+ sol -= (y3 + d3) * (x2 - y2 - d2) * G00(-0.5, x1, x2, d1, y2 + d2,\r
+ -y3 - d3);\r
+\r
+ if ((y3 + d3) != 0)\r
+ sol += (y3 + d3) * g0(0.5, x1, d1,\r
+ sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + (y3 + d3) * (y3 + d3)));\r
+\r
+ return sol / 2.;\r
+}\r
+\r
double quadInt(\r
double(*f)(double, double, double, double, double, double, double),\r
- double a, double b, double c, double d, double r, double t,\r
- double u, 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,v,d1,d2,d3)+f(b,d,r,u,d1,d2,d3)\r
- -f(b,c,t,v,d1,d2,d3)+f(b,c,t,u,d1,d2,d3)+f(b,c,r,v,d1,d2,d3)-f(b,c,r,u,d1,d2,d3)\r
- -f(a,d,t,v,d1,d2,d3)+f(a,d,t,u,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,v,d1,d2,d3)+f(a,c,r,u,d1,d2,d3);\r
-/*\r
- return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f(\r
- k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1,\r
- s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2,\r
- t1, l2, d1, d2, d3) - f(k1, s2, t1, t2, d1, d2, d3) - f(s1, k2, l1,\r
- l2, d1, d2, d3) + f(s1, k2, l1, t2, d1, d2, d3) + f(s1, k2, t1, l2,\r
- d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1,\r
- d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2,\r
- d3) + f(s1, s2, t1, t2, d1, d2, d3);\r
-\r
- */\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
+ /*\r
+ return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f(\r
+ k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1,\r
+ s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2,\r
+ t1, l2, d1, d2, d3) - f(k1, s2, t1, t2, d1, d2, d3) - f(s1, k2, l1,\r
+ l2, d1, d2, d3) + f(s1, k2, l1, t2, d1, d2, d3) + f(s1, k2, t1, l2,\r
+ d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1,\r
+ d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2,\r
+ d3) + f(s1, s2, t1, t2, d1, d2, d3);\r
+\r
+ */\r
}\r
\r
double quad0Int(\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)-f(b,d,t,0,d1,d2,d3)-f(b,d,0,v,d1,d2,d3)+f(b,d,0,0,d1,d2,d3)\r
- -f(b,0,t,v,d1,d2,d3)+f(b,0,t,0,d1,d2,d3)+f(b,0,0,v,d1,d2,d3)-f(b,0,0,0,d1,d2,d3)\r
- -f(0,d,t,v,d1,d2,d3)+f(0,d,t,0,d1,d2,d3)+f(0,d,0,v,d1,d2,d3)-f(0,d,0,0,d1,d2,d3)\r
- +f(0,0,t,v,d1,d2,d3)-f(0,0,t,0,d1,d2,d3)-f(0,0,0,v,d1,d2,d3)+f(0,0,0,0,d1,d2,d3);\r
+ return f(b, d, t, v, d1, d2, d3) - f(b, d, t, 0, d1, d2, d3) - f(b, d, 0,\r
+ v, d1, d2, d3) + f(b, d, 0, 0, d1, d2, d3) - f(b, 0, t, v, d1, d2,\r
+ d3) + f(b, 0, t, 0, d1, d2, d3) + f(b, 0, 0, v, d1, d2, d3) - f(b,\r
+ 0, 0, 0, d1, d2, d3) - f(0, d, t, v, d1, d2, d3) + f(0, d, t, 0,\r
+ d1, d2, d3) + f(0, d, 0, v, d1, d2, d3) - f(0, d, 0, 0, d1, d2, d3)\r
+ + f(0, 0, t, v, d1, d2, d3) - f(0, 0, t, 0, d1, d2, d3) - f(0, 0,\r
+ 0, v, d1, d2, d3) + f(0, 0, 0, 0, d1, d2, d3);\r
\r
}\r
\r
-\r
-\r
-\r
-\r
double slpADLO(double y1, double y2, double x1, double x2, double a) {\r
double G3 = 0;\r
double gL = 0;\r
return (y - x) * pow((x * x + y * y + a * a - 2 * x * y), -0.5)\r
/ (a * a);\r
default:\r
- printf("p must be either 0, -1/2, -1 or -3/2.");\r
+ printf("p must be either 0, -1/2, -1 or -3/2. (%.2f)\n", p);\r
return NAN;\r
}\r
}\r
x1, yd1, tmp3) - tmp2 * compute_g0(-0.5, x2, yd2, tmp4) + tmp5 / 3.;\r
return rval;\r
}\r
+\r
+double FLO_perp(double x1, double x2, double y2, double y3, double delta1,\r
+ double delta2, double a) {\r
+ double xd1 = x1 - delta1;\r
+ double xd2 = x2 - delta2;\r
+ double yd2 = y2 - delta2;\r
+ double yd3 = y3 - a;\r
+ double tmp2 = x2 - y2 - delta2;\r
+ double rval = 0;\r
+\r
+ rval = xd1 * compute_g0(y3, a, sqrt(xd1 * xd1 + tmp2 * tmp2), -0.5) + yd3\r
+ * compute_g0(x1, delta1, sqrt(tmp2 * tmp2 + yd3 * yd3), -0.5) - xd1\r
+ * tmp2 * slpADLO(x2, y3, yd2, a, xd1) - yd3 * tmp2 * slpADLO(x1,\r
+ x2, delta1, yd2, -y3 - a) - slpADLO(y3, x1, -a, delta1, tmp2);\r
+ return rval;\r
+}\r
double G00(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,y1,y2,d1,d2,d3);
-//double F_ort(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);
double slpADLO(double, double, double, double, double);
double compute_g0(double, double, double, double);
double FLO_plane(double, double, double, double, double, double, double);
+double FLO_perp(double, double, double, double, double, double, double);
// sol = quadInt((F_par/F_ort),a,b,c,d,r,t,u,v,d1,d2,d3);
double quadInt(
+++ /dev/null
-#include <iostream>\r
-#include <cmath>\r
-#include <cassert>\r
-#include <stdlib.h>\r
-\r
-#include "SLPrecangle.hpp"\r
-\r
-#define my_PI 3.141592653589793\r
-#define EPS 0.02\r
-\r
-using namespace std;\r
-\r
-int sign(double);\r
-double arsinh(double);\r
-\r
-// sol = quadInt((F_par/F_ort),s1,s2,k1,k2,t1,t2,l1,l2,d1,d2,d3);\r
-double quadInt(\r
- double(*)(double, double, double, double, double, double, double),\r
- double, double, double, double, double, double, double, double, double,\r
- double, double);\r
-\r
-void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
- //sicherheitsabfragen zu Datengroessen\r
- if (nrhs != 2)\r
- mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx3))");\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 != 3)\r
- mexErrMsgTxt("expected elements (Mx3)");\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 * 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
-\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
-\r
- double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
- double tmp;\r
-\r
- int rx, ry;\r
-\r
- //Ausrechnen\r
- for (int j = 0; j < em; ++j) {\r
- x1[0] = C[(int) E[j] - 1];\r
- x1[1] = C[cm + (int) E[j] - 1];\r
- x1[2] = C[2 * cm + (int) E[j] - 1];\r
-\r
- x2[0] = C[(int) E[em + j] - 1];\r
- x2[1] = C[cm + (int) E[em + j] - 1];\r
- x2[2] = C[2 * cm + (int) E[em + j] - 1];\r
-\r
- x3[0] = C[(int) E[2 * em + j] - 1];\r
- x3[1] = C[cm + (int) E[2 * em + j] - 1];\r
- x3[2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
-\r
- xn[0] = (x2[1] - x1[1]) * (x3[2] - x1[2]) - (x2[2] - x1[2]) * (x3[1]\r
- - x1[1]);\r
- xn[1] = (x2[2] - x1[2]) * (x3[0] - x1[0]) - (x2[0] - x1[0]) * (x3[2]\r
- - x1[2]);\r
- xn[2] = (x2[0] - x1[0]) * (x3[1] - x1[1]) - (x2[1] - x1[1]) * (x3[0]\r
- - x1[0]);\r
-\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
- // printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);\r
-\r
- for (int k = 0; k < em; ++k) {\r
- y1[0] = C[(int) E[k] - 1];\r
- y1[1] = C[cm + (int) E[k] - 1];\r
- y1[2] = C[2 * cm + (int) E[k] - 1];\r
-\r
- y2[0] = C[(int) E[em + k] - 1];\r
- y2[1] = C[cm + (int) E[em + k] - 1];\r
- y2[2] = C[2 * cm + (int) E[em + k] - 1];\r
-\r
- y3[0] = C[(int) E[2 * em + k] - 1];\r
- y3[1] = C[cm + (int) E[2 * em + k] - 1];\r
- y3[2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
-\r
- yn[0] = (y2[1] - y1[1]) * (y3[2] - y1[2]) - (y2[2] - y1[2])\r
- * (y3[1] - y1[1]);\r
- yn[1] = (y2[2] - y1[2]) * (y3[0] - y1[0]) - (y2[0] - y1[0])\r
- * (y3[2] - y1[2]);\r
- yn[2] = (y2[0] - y1[0]) * (y3[1] - y1[1]) - (y2[1] - y1[1])\r
- * (y3[0] - y1[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
- d[0] = y1[0] - x1[0];\r
- d[1] = y1[1] - x1[1];\r
- d[2] = y1[2] - x1[2];\r
-\r
- if (rx = ry) {\r
- if (rx = 2) {\r
- //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x1[0],x1[1],x2[0],x3[1],y1[0],y1[1],y2[0],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 = quadInt(FLO_plane, x1[0], x1[1], x2[0], x3[1], y1[0],\r
- y1[1], y2[0], y3[1], d[0], d[1], d[2]);\r
- // printf("%d%d|%.2f\n",j,k,tmp);\r
- A[(k * em) + j] = 1. / (4 * my_PI) * tmp;\r
-\r
- } else\r
- A[(k * em) + j] = 0;\r
- } else {\r
- A[(k * em) + j] = 0;\r
- }\r
-\r
- // Vorbereiten der DATEN\r
- // ej = coordinates(elements(j,:)',:);\r
- // ek = coordinates(elements(k,:)',:);\r
- // d = ek(1,:) - ej(1,:);\r
-\r
- // AUSRECHNEN\r
- // A[(k*em)+j] = 1./(4*my_PI);// *\r
- // quadInt(F_par,\r
- // ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));\r
-\r
- }\r
-\r
- }\r
- //Rueckgabe (eventuell zurueck schreiben)\r
- // mxFree(x1);\r
- //mxFree(x2);\r
- //mxFree(x3);\r
- //mxFree(y1);\r
- //mxFree(y2);\r
- //mxFree(y3);\r
- //mxFree(xn);\r
- //mxFree(yn);\r
- //mxFree(d);\r
-\r
- return;\r
-}\r
-\r
-int inline sign(double x) {\r
- return x > 0 ? 1 : (x < 0 ? -1 : 0);\r
-}\r
-\r
-double inline arsinh(double x) {\r
- return log(x + sqrt(x * x + 1));\r
-}\r
-\r
-//y-x muss != 0 sein\r
-double g0(double p, double y, double x, double l) {\r
- //printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);\r
-\r
- double sol = 0;\r
-\r
- if (l != 0) {\r
- if (p == 0.5) {\r
- sol = (y - x) / 2 * sqrt((y - x) * (y - x) + l * l) + l * l / 2\r
- * arsinh((y - x) / fabs(l));\r
- // printf("%.2f |",sol);\r
- } else if (p == 0)\r
- sol = y - x;\r
- else if (p == -0.5)\r
- sol = asinh((y - x) / fabs(l));\r
- else if (p == -1)\r
- sol = atan((y - x) / fabs(l));\r
- else if (p == -1.5)\r
- sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l));\r
- else\r
- sol = (y - x) * pow((y - x) * (y - x) + l * l, p) + 2 * p * l * l\r
- * g0(p - 1, y, x, l) / (2 * p + 1);\r
- } else {\r
- if (p == -0.5)\r
- sol = sign(y - x) * log(fabs(y - x));\r
- else\r
- sol = (y - x) * pow(fabs(y - x), 2 * p) / (2 * p + 1);\r
- }\r
-\r
- return sol;\r
-}\r
-\r
-double G00(double p, double y1, double y2, double x1, double x2, double l) {\r
- // printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);\r
-\r
- double sol = 0;\r
- if (p == -1.5) {\r
- if (l == 0) {\r
- sol = -sqrt((y1 - x1) * (y1 - x1) + (y2 - x2) * (y2 - x2)) / ((y1\r
- - x1) * (y2 - x2));\r
- } else {\r
- sol = sign((y1 - x1) * (y2 - x2)) / (2 * fabs(l)) * acos(\r
- -2 * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2) / (((y1\r
- - x1) * (y1 - x1) + l * l) * ((y2 - x2) * (y2 - x2)\r
- + l * l)) + 1);\r
- }\r
-\r
- } else if (p == -0.5) {\r
- if (l != 0)\r
- sol = 2 * p * l * l * G00(p - 1, y1, y2, x1, x2, l);\r
- if ((y1 - x1) != 0)\r
- sol += (y1 - x1) * g0(p, y2, x2,\r
- sqrt((y1 - x1) * (y1 - x1) + l * l));\r
- if ((y2 - x2) != 0)\r
- sol += (y2 - x2) * g0(p, y1, x1,\r
- sqrt((y2 - x2) * (y2 - x2) + l * l));\r
- sol /= 2 * p + 2;\r
- } else {\r
- mexErrMsgTxt("no case for p defined");\r
- }\r
-\r
- return sol;\r
-}\r
-\r
-double F_par(double x1, double x2, double y1, double y2, double d1, double d2,\r
- double d3) {\r
-\r
- // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
- double sol = (x1 - y1 - d1) * (x2 - y2 - d2);\r
-\r
- if (sol != 0)\r
- sol *= slpADLO(x1, x2, y1 + d1, y2 + d2, d3);\r
-\r
- if ((x1 - y1 - d1) != 0)\r
- sol -= (x1 - y1 - d1) * compute_g0(0.5, x1, y1 + d1,\r
- sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));\r
-\r
- if ((x2 - y2 - d2) != 0)\r
- sol -= (x2 - y2 - d2) * compute_g0(0.5, x2, y2 + d2,\r
- sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));\r
-\r
- double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
- - d2) + d3 * d3);\r
- sol += 1. / 3 * hlp * sqrt(hlp);\r
- return sol;\r
-}\r
-\r
-double inline quadInt(\r
- double(*f)(double, double, double, double, double, double, double),\r
- double s1, double s2, double k1, double k2, double t1, double t2,\r
- double l1, double l2, double d1, double d2, double d3) {\r
-\r
- return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f(\r
- k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1,\r
- s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2,\r
- t1, l2, d1, d2, d3) - f(k1, s2, t1, t2, d1, d2, d3) - f(s1, k2, l1,\r
- l2, d1, d2, d3) + f(s1, k2, l1, t2, d1, d2, d3) + f(s1, k2, t1, l2,\r
- d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1,\r
- d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2,\r
- d3) + f(s1, s2, t1, t2, d1, d2, d3);\r
-}\r
-\r
-double slpADLO(double y1, double y2, double x1, double x2, double a) {\r
- double G3 = 0;\r
- double gL = 0;\r
- double gK = 0;\r
- double tmp;\r
-\r
- tmp = y1 - x1;\r
- if (fabs(tmp) >= EPS * y1) {\r
- tmp = sqrt(y1 * y1 + x1 * x1 + a * a - 2 * y1 * x1);\r
- gL = compute_g0(-0.5, y2, x2, tmp);\r
- }\r
- tmp = y2 - x2;\r
- if (fabs(tmp) >= EPS * y2) {\r
- tmp = sqrt(y2 * y2 + x2 * x2 + a * a - 2 * y2 * x2);\r
- gK = compute_g0(-0.5, y1, x1, tmp);\r
- }\r
- if (fabs(a * a) > EPS) {\r
- if ((y1 - x1) * (y2 - x2) * a >= 0)\r
- tmp = 1.;\r
- else\r
- tmp = -1.;\r
-\r
- G3 = tmp * acos(\r
- (-2. * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2)) / (((y1\r
- - x1) * (y1 - x1) + a * a) * ((y2 - x2) * (y2 - x2) + a\r
- * a)) + 1.) / (2. * a);\r
- }\r
-\r
- return (y1 - x1) * gL + (y2 - x2) * gK - a * a * G3;\r
-}\r
-\r
-double compute_g0(double p, double y, double x, double a) {\r
- int sp = (int) 2 * (p - EPS); // MK (p-EPS) instead of (p+EPS)\r
- //printf("\n compute_g0, p = %lf, sp = %d\n",p,sp);\r
- assert(\r
- p == 0 || (p == -0.5) || (p == -1) || (p == -1.5) || (fabs(a)\r
- <= EPS));\r
- if (fabs(a) <= EPS) {\r
- // printf("\n a < eps\n");\r
- switch (sp) {\r
- case 0:\r
- return y - x;\r
- case -1:\r
- return log(fabs(x - y)) * (y - x) / fabs(y - x);\r
- case -2:\r
- return -(y - x) / fabs(y - x) / fabs(y - x);\r
- case -3:\r
- return -0.5 * (y - x) / fabs(y - x) / fabs(y - x) / fabs(y - x);\r
- }\r
- } else {\r
- // printf("\n a > eps\n");\r
- switch (sp) {\r
- case 0:\r
- return y - x;\r
- case -1:\r
- return asinh((y - x) / fabs(a));\r
- case -2:\r
- return atan((y - x) / fabs(a));\r
- case -3:\r
- return (y - x) * pow((x * x + y * y + a * a - 2 * x * y), -0.5)\r
- / (a * a);\r
- default:\r
- printf("p must be either 0, -1/2, -1 or -3/2.");\r
- return NAN;\r
- }\r
- }\r
-}\r
-\r
-double FLO_plane(double x1, double x2, double y1, double y2, double delta1,\r
- double delta2, double a) {\r
- double yd1 = y1 + delta1;\r
- double yd2 = y2 + delta2;\r
- double tmp1 = x1 - y1 - delta1;\r
- double tmp2 = x2 - y2 - delta2;\r
- double tmp3 = sqrt(tmp2 * tmp2 + a * a);\r
- double tmp4 = sqrt(tmp1 * tmp1 + a * a);\r
- double tmp5 = pow(tmp1 * tmp1 + tmp2 * tmp2 + a * a, 3. / 2.);\r
- double rval = 0;\r
-\r
- rval = tmp1 * tmp2 * slpADLO(x1, x2, yd1, yd2, a) - tmp1 * compute_g0(-0.5,\r
- x1, yd1, tmp3) - tmp2 * compute_g0(-0.5, x2, yd2, tmp4) + tmp5 / 3.;\r
- return rval;\r
-}\r
+++ /dev/null
-#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 g0(double, double, double, double);
-// sol = G00(p,y1,y2,x1,x2,l);
-double G00(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,y1,y2,d1,d2,d3);
-//double F_ort(double, double, double, double, double, double, double);
-
-
-double slpADLO(double, double, double, double, double);
-double compute_g0(double, double, double, double);
-double FLO_plane(double, double, double, double, double, double, double);
-
-#endif
--- /dev/null
+#include <iostream>\r
+#include <cmath>\r
+#include <cassert>\r
+#include <stdlib.h>\r
+\r
+#include "SLPrectangle.hpp"\r
+\r
+#define my_PI 3.141592653589793\r
+#define EPS 0.00001\r
+\r
+using namespace std;\r
+\r
+int sign(double);\r
+//double arsinh(double);\r
+\r
+\r
+int inline sign(double x) {\r
+ return x > 0 ? 1 : (x < 0 ? -1 : 0);\r
+}\r
+/*\r
+ double inline arsinh(double x) {\r
+ return log(x + sqrt(x * x + 1));\r
+ }\r
+ */\r
+\r
+//y-x muss != 0 sein\r
+double g0(double p, double y, double x, double l) {\r
+ //printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);\r
+\r
+ double sol = 0;\r
+\r
+ if (l != 0) {\r
+ if (p == 0.5) {\r
+ sol = (y - x) / 2 * sqrt((y - x) * (y - x) + l * l) + l * l / 2\r
+ * asinh((y - x) / fabs(l));\r
+ // printf("%.2f |",sol);\r
+ } else if (p == 0)\r
+ sol = y - x;\r
+ else if (p == -0.5)\r
+ sol = asinh((y - x) / fabs(l));\r
+ else if (p == -1)\r
+ sol = atan((y - x) / fabs(l)) / l;\r
+ else if (p == -1.5)\r
+ sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l));\r
+ else\r
+ sol = (y - x) * pow((y - x) * (y - x) + l * l, p) + 2 * p * l * l\r
+ * g0(p - 1, y, x, l) / (2 * p + 1);\r
+ } else {\r
+ if (p == -0.5)\r
+ sol = sign(y - x) * log(fabs(y - x));\r
+ else\r
+ sol = (y - x) * pow(fabs(y - x), 2 * p) / (2 * p + 1);\r
+ }\r
+\r
+ return sol;\r
+}\r
+\r
+double G00(double p, double y1, double y2, double x1, double x2, double l) {\r
+ // printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);\r
+\r
+ double sol = 0;\r
+ if (p == -1.5) {\r
+ if (l == 0) {\r
+ sol = -sqrt((y1 - x1) * (y1 - x1) + (y2 - x2) * (y2 - x2)) / ((y1\r
+ - x1) * (y2 - x2));\r
+ } else {\r
+ sol = sign((y1 - x1) * (y2 - x2)) / (2 * fabs(l)) * acos(\r
+ -2 * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2) / (((y1\r
+ - x1) * (y1 - x1) + l * l) * ((y2 - x2) * (y2 - x2)\r
+ + l * l)) + 1);\r
+ }\r
+\r
+ } else if ((pt > 0) && ((int) pt == pt)) {\r
+ if (l != 0)\r
+ sol = 2 * p * l * l * G00(p - 1, y1, y2, x1, x2, l);\r
+ if ((y1 - x1) != 0)\r
+ sol += (y1 - x1) * g0(p, y2, x2,\r
+ sqrt((y1 - x1) * (y1 - x1) + l * l));\r
+ if ((y2 - x2) != 0)\r
+ sol += (y2 - x2) * g0(p, y1, x1,\r
+ sqrt((y2 - x2) * (y2 - x2) + l * l));\r
+ sol /= 2 * p + 2;\r
+ } else {\r
+ sol = NAN;\r
+ printf("warning in G00: no case for p=%.2f defined\n", p);\r
+ }\r
+\r
+ return sol;\r
+}\r
+\r
+double F_par(double x1, double x2, double y1, double y2, double d1, double d2,\r
+ double d3) {\r
+\r
+ // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
+ double sol = (x1 - y1 - d1) * (x2 - y2 - d2);\r
+\r
+ if (sol != 0)\r
+ sol *= G00(-0.5, x1, x2, y1 + d1, y2 + d2, d3);\r
+\r
+ if ((x1 - y1 - d1) != 0)\r
+ sol -= (x1 - y1 - d1) * g0(0.5, x1, y1 + d1,\r
+ sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));\r
+\r
+ if ((x2 - y2 - d2) != 0)\r
+ sol -= (x2 - y2 - d2) * g0(0.5, x2, y2 + d2,\r
+ sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));\r
+\r
+ double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
+ - d2) + d3 * d3);\r
+ sol += 1. / 3 * hlp * sqrt(hlp);\r
+ return sol;\r
+}\r
+\r
+double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,\r
+ double d3) {\r
+\r
+ // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y2,y3,d1,d2,d3);\r
+ double sol = -G00(0.5, y3, x1, -d3, d1, x2 - y2 - d2);\r
+\r
+ if ((x1 - d1) * (x2 - y2 - d2) != 0)\r
+ sol -= (x1 - d1) * (x2 - y2 - d2) * G00(-0.5, x2, y3, y2 + d2, -d3,\r
+ x1 - d1);\r
+\r
+ if ((x1 - d1) != 0)\r
+ sol += (x1 - d1) * g0(0.5, y3, -d3,\r
+ sqrt((x1 - d1) * (x1 - d1) + (x2 - y2 - d2) * (x2 - y2 - d2)));\r
+\r
+ if ((y3 + d3) * (x2 - y2 - d2) != 0)\r
+ sol -= (y3 + d3) * (x2 - y2 - d2) * G00(-0.5, x1, x2, d1, y2 + d2,\r
+ -y3 - d3);\r
+\r
+ if ((y3 + d3) != 0)\r
+ sol += (y3 + d3) * g0(0.5, x1, d1,\r
+ sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + (y3 + d3) * (y3 + d3)));\r
+\r
+ return sol / 2.;\r
+}\r
+\r
+double quadInt(\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
+}\r
+\r
+double quad0Int(\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) - f(b, d, t, 0, d1, d2, d3) - f(b, d, 0,\r
+ v, d1, d2, d3) + f(b, d, 0, 0, d1, d2, d3) - f(b, 0, t, v, d1, d2,\r
+ d3) + f(b, 0, t, 0, d1, d2, d3) + f(b, 0, 0, v, d1, d2, d3) - f(b,\r
+ 0, 0, 0, d1, d2, d3) - f(0, d, t, v, d1, d2, d3) + f(0, d, t, 0,\r
+ d1, d2, d3) + f(0, d, 0, v, d1, d2, d3) - f(0, d, 0, 0, d1, d2, d3)\r
+ + f(0, 0, t, v, d1, d2, d3) - f(0, 0, t, 0, d1, d2, d3) - f(0, 0,\r
+ 0, v, d1, d2, d3) + f(0, 0, 0, 0, d1, d2, d3);\r
+\r
+}\r
+\r
--- /dev/null
+#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 g0(double, double, double, double);
+// sol = G00(p,y1,y2,x1,x2,l);
+double G00(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 = quadInt((F_par/F_ort),a,b,c,d,r,t,u,v,d1,d2,d3);
+double quadInt(
+ double(*)(double, double, double, double, double, double, double),
+ double, double, double, double, double, double, double, double, double,
+ double, double);
+
+// sol = quad0Int((F_par/F_ort),b,d,t,v,d1,d2,d3);
+double quad0Int(
+ double(*)(double, double, double, double, double, double, double),
+ double, double, double, double, double, double, double);
+
+
+#endif
#include "../spaces/P0Space.hpp"
#include "TriangleIntegrator.hpp"
-#include "SLPrecangle.hpp"
-
+#include "SLPrectangle.hpp" // Bibliothek fuer achsenorientierte Rechtecke
#include <fstream>
using namespace laplace;
-SimpleLayerPotential::SimpleLayerPotential(const space::P0Space& space)
- : HLibOperator< space::P0Space, space::P0Space, SymmetricMatrix >(
- space, space, SimpleLayerPotential::computeEntry)
-{
+SimpleLayerPotential::SimpleLayerPotential(const space::P0Space& space) :
+ HLibOperator<space::P0Space, space::P0Space, SymmetricMatrix> (
+ space, space, SimpleLayerPotential::computeEntry) {
}
-SimpleLayerPotential::~SimpleLayerPotential()
-{
+SimpleLayerPotential::~SimpleLayerPotential() {
}
-template < class TMesh >
-double ComputeSimpleLayerPotentialEntry< TMesh >::exec(
- const typename TMesh::Prop::Triangle& f1,
- const typename TMesh::Prop::Triangle& f2) const
-{
- TriangleIntegrator integrator = TriangleIntegrator();
-
- if ( f1.computeArea() < f2.computeArea() )
- {
- return integrator.integrate( f1, f2,
- computeInnerIntegralForTriangles);
- }
- else
- {
- return integrator.integrate( f2, f1,
- computeInnerIntegralForTriangles);
- }
+template<class TMesh>
+double ComputeSimpleLayerPotentialEntry<TMesh>::exec(
+ const typename TMesh::Prop::Triangle& f1,
+ const typename TMesh::Prop::Triangle& f2) const {
+ TriangleIntegrator integrator = TriangleIntegrator();
+
+ if (f1.computeArea() < f2.computeArea()) {
+ return integrator.integrate(f1, f2, computeInnerIntegralForTriangles);
+ } else {
+ return integrator.integrate(f2, f1, computeInnerIntegralForTriangles);
+ }
}
-template < class TMesh >
-double ComputeSimpleLayerPotentialEntry< TMesh >::exec(
- const typename TMesh::Prop::Parallelogram& f1,
- const typename TMesh::Prop::Parallelogram& f2) const
-{
-
- const typename TMesh::Vertex& a1 = f1.getA();
- const typename TMesh::Vertex& a2 = f2.getA();
-
-
-
-
- return NAN;
+template<class TMesh>
+double ComputeSimpleLayerPotentialEntry<TMesh>::exec(
+ const typename TMesh::Prop::Parallelogram& f1,
+ const typename TMesh::Prop::Parallelogram& f2) const {
+
+ // Zwischenvariable fuer die Loesung (kann fuer geschwindigkeit entfernt werden)
+ double sol = NAN;
+
+ // Flaechenstueck 1
+ const typename TMesh::Vertex& xa = f1.getA() - f1.getB();
+ const typename TMesh::Vertex& xb = f1.getA() - f1.getD();
+ const typename TMesh::Vertex& xn = f1.getA() * f1.getD(); //todo Kreuzprodukt!?
+
+ // Flaechenstueck 2
+ const typename TMesh::Vertex& ya = f2.getA() - f1.getB();
+ const typename TMesh::Vertex& yb = f2.getA() - f1.getD();
+ const typename TMesh::Vertex& yn = f2.getA() * f2.getD(); //todo Kreuzprodukt!?
+
+ // Abstand zwischen den Stuecken (Achtung) abstand von den beiden "kleinsten" Ecken
+ const typename TMesh::Vertex& d = NULL; //todo new???
+
+ //Variablen fuer LageInformationen
+ int rx, rxa, rxb, ry, rya, ryb;
+
+ // Lage des ersten Elements ermitteln
+ if (xn[2] != 0)
+ rx = 2;
+ else if (xn[1] != 0)
+ rx = 1;
+ else
+ rx = 0;
+
+ if (xa[2] != 0)
+ rxa = 2;
+ else if (xa[1] != 0)
+ rxa = 1;
+ else
+ rxa = 0;
+
+ if (xb[2] != 0)
+ rxb = 2;
+ else if (xb[1] != 0)
+ rxb = 1;
+ else
+ rxb = 0;
+
+ //kleinste Ecke finden und fuer \delta verwenden
+ if (xa[rxa] > 0) {
+ if (xb[rxb] > 0) {
+ d = -f1.getA;
+ } else {
+ d = -f1.getD;
+ }
+ } else {
+ if (xb[rxb] > 0) {
+ d = -f1.getB;
+ } else {
+ d = -f1.getC;
+ }
+ }
+
+ // Lage des zweiten Elements ermitteln
+ if (yn[2] != 0)
+ ry = 2;
+ else if (yn[1] != 0)
+ ry = 1;
+ else
+ ry = 0;
+
+ if (ya[2] != 0)
+ rya = 2;
+ else if (ya[1] != 0)
+ rya = 1;
+ else
+ rya = 0;
+
+ if (yb[2] != 0)
+ ryb = 2;
+ else if (yb[1] != 0)
+ ryb = 1;
+ else
+ ryb = 0;
+
+ //kleinste Ecke finden und fuer \delta verwenden
+ if (ya[rya] > 0) {
+ if (yb[ryb] > 0) {
+ d += f2.getA; //todo += definiert??? sonst d = d + ...
+ } else {
+ d += f2.getD; //todo += definiert??? sonst d = d + ...
+ }
+ } else {
+ if (yb[ryb] > 0) {
+ d += f2.getB; //todo += definiert??? sonst d = d + ...
+ } else {
+ d += f2.getC; //todo += definiert??? sonst d = d + ...
+ }
+ }
+
+ if (rx == ry) { // Flaechen sind parallel
+ if (rxa == rya) { // Elemente zeigen in gleiche Richtung
+ sol = quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]),
+ fabs(yb[rxb]), d[rxa], d[rxb], d[rx]);
+ } else { // Elemente zeigen nicht in gleiche Richtung
+ sol = quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[rxa]),
+ fabs(ya[rxb]), d[rxa], d[rxb], d[rx]);
+ }
+ } else { // Flaechen sind orthogonal
+ if (rxa == rya) { //Elmente richtig zuordnen
+ sol = quad0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]), fabs(ya[rya]),
+ fabs(yb[ryb]), d[rxb], d[rxa], d[rx]);
+ } else if (rxa == ryb) {
+ sol = quad0Int(F_ort, fabs(xb[rxb]), fabs(xa[rxa]), fabs(yb[ryb]),
+ fabs(ya[rya]), d[rxb], d[rxa], d[rx]);
+ } else if (rxb == rya) {
+ sol = quad0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rya]),
+ fabs(yb[ryb]), d[rxa], d[rxb], d[rx]);
+ } else {
+ sol = quad0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[ryb]),
+ fabs(ya[rya]), d[rxa], d[rxb], d[rx]);
+ }
+ sol /= 2.;
+ }
+
+ return sol / (4. * PI);
}
-double SimpleLayerPotential::computeEntry(
- const space::P0Space::Element &e1,
- const space::P0Space::Element &e2)
-{
- typedef boundary_mesh::DoubleFaceDispatcher<
- ComputeSimpleLayerPotentialEntry< BoundaryMesh >, BoundaryMesh, double >
- ComputeEntry;
+double SimpleLayerPotential::computeEntry(const space::P0Space::Element &e1,
+ const space::P0Space::Element &e2) {
+ typedef boundary_mesh::DoubleFaceDispatcher<
+ ComputeSimpleLayerPotentialEntry<BoundaryMesh> , BoundaryMesh,
+ double> ComputeEntry;
- return ComputeEntry::Exec( e1.getSupport(), e2.getSupport(),
- ComputeSimpleLayerPotentialEntry< BoundaryMesh >() );
+ return ComputeEntry::Exec(e1.getSupport(), e2.getSupport(),
+ ComputeSimpleLayerPotentialEntry<BoundaryMesh> ());
}
-template < class TMesh >
-double
-ComputeSimpleLayerPotentialEntry< TMesh >::computeInnerIntegralForTriangles(
- const typename TMesh::Prop::Triangle& tau,
- const typename TMesh::UVertex& x)
-{
- BoundaryMesh::Face::ConstVertexIterator vIt = tau.begin_v();
- const BoundaryMesh::UVertex& x1 = *vIt; ++vIt;
- const BoundaryMesh::UVertex& x2 = *vIt; ++vIt;
- const BoundaryMesh::UVertex& x3 = *vIt; ++vIt;
-
- // determine a local orthogonal coordinate system (r1,r2,n) on the triangle tau
- double t_tau = (x3-x2).norm();
- BoundaryMesh::UVertex r2 = (x3-x2)/t_tau;
- double t_star=(x1-x2)*r2;
-
- BoundaryMesh::UVertex r1 = (x2+t_star*r2-x1);
- double s_tau = r1.norm();
- r1 = r1 / s_tau;
-
- BoundaryMesh::UVertex n = BoundaryMesh::Vertex::crossProduct(r1,r2);
-
- // determine necessary information
- double sx = (x-x1)*r1;
- double tx = (x-x1)*r2;
- double ux = (x-x1)*n;
- double a1 = -t_star/s_tau;
- double a2 = (t_tau-t_star)/s_tau;
-
- // analytic computation of inner integral
- return 0.25*(F(s_tau,a2,sx,tx,ux) - F(0,a2,sx,tx,ux) - F(s_tau,a1,sx,tx,ux) + F(0,a1,sx,tx,ux))/M_PI;
+template<class TMesh>
+double ComputeSimpleLayerPotentialEntry<TMesh>::computeInnerIntegralForTriangles(
+ const typename TMesh::Prop::Triangle& tau,
+ const typename TMesh::UVertex& x) {
+ BoundaryMesh::Face::ConstVertexIterator vIt = tau.begin_v();
+ const BoundaryMesh::UVertex& x1 = *vIt;
+ ++vIt;
+ const BoundaryMesh::UVertex& x2 = *vIt;
+ ++vIt;
+ const BoundaryMesh::UVertex& x3 = *vIt;
+ ++vIt;
+
+ // determine a local orthogonal coordinate system (r1,r2,n) on the triangle tau
+ double t_tau = (x3 - x2).norm();
+ BoundaryMesh::UVertex r2 = (x3 - x2) / t_tau;
+ double t_star = (x1 - x2) * r2;
+
+ BoundaryMesh::UVertex r1 = (x2 + t_star * r2 - x1);
+ double s_tau = r1.norm();
+ r1 = r1 / s_tau;
+
+ BoundaryMesh::UVertex n = BoundaryMesh::Vertex::crossProduct(r1, r2);
+
+ // determine necessary information
+ double sx = (x - x1) * r1;
+ double tx = (x - x1) * r2;
+ double ux = (x - x1) * n;
+ double a1 = -t_star / s_tau;
+ double a2 = (t_tau - t_star) / s_tau;
+
+ // analytic computation of inner integral
+ return 0.25 * (F(s_tau, a2, sx, tx, ux) - F(0, a2, sx, tx, ux) - F(s_tau,
+ a1, sx, tx, ux) + F(0, a1, sx, tx, ux)) / M_PI;
}
-template < class TMesh >
-double
-ComputeSimpleLayerPotentialEntry< TMesh >::F(
- double s, double a, double sx, double tx, double ux)
-{
- double p = (a*tx+sx)/(1+a*a);
- double q = ux*ux + (tx-a*sx)*(tx-a*sx)/(1+a*a);
-
- return (s-sx)*log(a*s - tx + sqrt((s-sx)*(s-sx) + (a*s-tx)*(a*s-tx) + ux*ux)) - s
- + (a*sx-tx)*log(sqrt(1+a*a)*(s-p)+sqrt((1+a*a)*(s-p)*(s-p) + q*q))/sqrt(q+a*a)
- + 2*ux*atan(((q-(a*sx-tx)/(1+a*a))*sqrt((1+a*a)*(s-p)*(s-p)+q*q) + (a*s-tx-q)*q)/((s-p)*ux));
+template<class TMesh>
+double ComputeSimpleLayerPotentialEntry<TMesh>::F(double s, double a,
+ double sx, double tx, double ux) {
+ double p = (a * tx + sx) / (1 + a * a);
+ double q = ux * ux + (tx - a * sx) * (tx - a * sx) / (1 + a * a);
+
+ return (s - sx) * log(
+ a * s - tx
+ + sqrt(
+ (s - sx) * (s - sx) + (a * s - tx) * (a * s - tx)
+ + ux * ux)) - s + (a * sx - tx) * log(
+ sqrt(1 + a * a) * (s - p) + sqrt(
+ (1 + a * a) * (s - p) * (s - p) + q * q)) / sqrt(q + a * a)
+ + 2 * ux * atan(
+ ((q - (a * sx - tx) / (1 + a * a)) * sqrt(
+ (1 + a * a) * (s - p) * (s - p) + q * q) + (a * s
+ - tx - q) * q) / ((s - p) * ux));
}
-void fill_supermatrix_aca_OLD(supermatrix* A, const cluster* row, const cluster* col, void *data,
- double (*compute_entry)(int i, int j, void* data), double aca_eps)
-{
-// std::cerr << "*** fill_supermatrix_aca" << std::endl;
-
- fullmatrix* f = A->f;
- rkmatrix* r = A->r;
- int i, j;
- int rank;
-
- if ( f!= NULL ) { // fill fullmatrix block
-// std::cerr << "*** full matrix block" << std::endl;
- for ( i=0; i<row->size; i++ ){
- for ( j=0; j<col->size; j++ )
- f->e[i+j*row->size] = compute_entry(row->start+i,col->start+j, data);
- }
-
- return;
- }
-
- if ( r!= NULL ) { // fill rk-block with aca
-// std::cerr << "*** rank k matrix block" << std::endl;
- if ( r->rows > r-> cols )
- rank = r->cols;
- else
- rank = r->rows;
-
- if ( rank > ACA_KMAX )
- rank = ACA_KMAX;
-
- r->k = rank;
- r->a = allocate_matrix(r->rows,rank);
- r->b = allocate_matrix(r->cols,rank);
-
- r->kt = newaca_fill_block(r->a, r->b, r->rows, r->cols,
- row->start, col->start,
- compute_entry, data,
- rank, aca_eps, ACA_STRATEGY);
-
- r->a = (double*) realloc(r->a,r->rows*r->kt*sizeof(double));
- r->b = (double*) realloc(r->b,r->cols*r->kt*sizeof(double));
- r->k = r->kt;
-
- return;
- }
-
- if(row->sons>0) {
- if(col->sons>0) {
- for(int i=0;i<row->sons;++i)
- for(int j=0;j<col->sons;++j) {
- int s=i+j*row->sons;
- fill_supermatrix_aca_OLD(A->s[s],row->son[i],col->son[j],data,compute_entry,aca_eps);
- }
- } else {
- for(int i=0;i<row->sons;++i)
- fill_supermatrix_aca_OLD(A->s[i],row->son[i],col,data,compute_entry,aca_eps);
- }
- } else {
- for(int i=0;i<col->sons;++i)
- fill_supermatrix_aca_OLD(A->s[i],row,col->son[i],data,compute_entry,aca_eps);
- }
+void fill_supermatrix_aca_OLD(supermatrix* A, const cluster* row,
+ const cluster* col, void *data,
+ double(*compute_entry)(int i, int j, void* data), double aca_eps) {
+ // std::cerr << "*** fill_supermatrix_aca" << std::endl;
+
+ fullmatrix* f = A->f;
+ rkmatrix* r = A->r;
+ int i, j;
+ int rank;
+
+ if (f != NULL) { // fill fullmatrix block
+ // std::cerr << "*** full matrix block" << std::endl;
+ for (i = 0; i < row->size; i++) {
+ for (j = 0; j < col->size; j++)
+ f->e[i + j * row->size] = compute_entry(row->start + i,
+ col->start + j, data);
+ }
+
+ return;
+ }
+
+ if (r != NULL) { // fill rk-block with aca
+ // std::cerr << "*** rank k matrix block" << std::endl;
+ if (r->rows > r-> cols)
+ rank = r->cols;
+ else
+ rank = r->rows;
+
+ if (rank > ACA_KMAX)
+ rank = ACA_KMAX;
+
+ r->k = rank;
+ r->a = allocate_matrix(r->rows, rank);
+ r->b = allocate_matrix(r->cols, rank);
+
+ r->kt = newaca_fill_block(r->a, r->b, r->rows, r->cols, row->start,
+ col->start, compute_entry, data, rank, aca_eps, ACA_STRATEGY);
+
+ r->a = (double*) realloc(r->a, r->rows * r->kt * sizeof(double));
+ r->b = (double*) realloc(r->b, r->cols * r->kt * sizeof(double));
+ r->k = r->kt;
+
+ return;
+ }
+
+ if (row->sons > 0) {
+ if (col->sons > 0) {
+ for (int i = 0; i < row->sons; ++i)
+ for (int j = 0; j < col->sons; ++j) {
+ int s = i + j * row->sons;
+ fill_supermatrix_aca_OLD(A->s[s], row->son[i], col->son[j],
+ data, compute_entry, aca_eps);
+ }
+ } else {
+ for (int i = 0; i < row->sons; ++i)
+ fill_supermatrix_aca_OLD(A->s[i], row->son[i], col, data,
+ compute_entry, aca_eps);
+ }
+ } else {
+ for (int i = 0; i < col->sons; ++i)
+ fill_supermatrix_aca_OLD(A->s[i], row, col->son[i], data,
+ compute_entry, aca_eps);
+ }
}
\r
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
\r
- int i,j,k;\r
+ int i, j, k;\r
\r
//sicherheitsabfragen zu Datengroessen\r
if (nrhs != 2)\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
+ for (i = 0; i < 3; ++i)\r
xa[i] = x1[i] - x0[i];\r
\r
- for (i = 0; i<3;++i)\r
+ for (i = 0; i < 3; ++i)\r
xb[i] = x3[i] - x0[i];\r
\r
//Kreuzprodukt\r
xn[1] = (xa[2]) * (xb[0]) - (xa[0]) * (xb[2]);\r
xn[2] = (xa[0]) * (xb[1]) - (xa[1]) * (xb[0]);\r
\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
\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
+ 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
+ } 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
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
+ for (i = 0; i < 3; ++i)\r
ya[i] = y1[i] - y0[i];\r
\r
- for (i = 0; i<3;++i)\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
\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
+ 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
+ } 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
+ } 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
+ } 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
+ // 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 %.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 = quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),\r
- fabs(ya[rxa]), fabs(yb[rxb]), d[rxa], d[rxb], d[rx]);\r
+ tmp\r
+ = quad0Int(F_par, 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
+ } 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 = quad0Int(F_par, fabs(xa[rxa]), fabs(xb[rxb]),\r
- fabs(yb[rxa]), fabs(ya[rxb]), d[rxa], d[rxb], d[rx]);\r
+ tmp\r
+ = quad0Int(F_par, 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
- A[(k * em) + j] = 1. / (4 * my_PI) * tmp;\r
\r
} else {\r
- A[(k * em) + j] = NAN;\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
+ = quad0Int(F_ort, 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
+ = quad0Int(F_ort, 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
+ = quad0Int(F_ort, 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
+ = quad0Int(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),\r
+ fabs(yb[ryb]), fabs(ya[rya]), d[rxa],\r
+ d[rxb], d[rx]);\r
+ }\r
+ A[(k * em) + j] = tmp;\r
+\r
}\r
+ A[(k * em) + j] = 1. / (4 * my_PI) * tmp;\r
\r
// Vorbereiten der DATEN\r
// ej = coordinates(elements(j,:)',:);\r
end
end
if(size(ele,2)==4)
-% disp 'schreibe in Datei'
-% dat = fopen(file,'w');
-%
-% for i=1:size(coo,1);
-% fprintf(dat,'v %.4f %.4f %.4f\n',coo(i,:));
-% end
-% fprintf(dat,'\n');
-% for i=1:size(ele,1);
-% fprintf(dat,'f %.0f %.0f %.0f %.0f\n',ele(i,:));
-% end
-%
-% fclose(dat);
+ disp 'schreibe in Datei'
+ dat = fopen(file,'w');
+
+ for i=1:size(coo,1);
+ fprintf(dat,'v %.4f %.4f %.4f\n',coo(i,:));
+ end
+ fprintf(dat,'\n');
+ for i=1:size(ele,1);
+ fprintf(dat,'f %.0f %.0f %.0f %.0f\n',ele(i,:));
+ end
+
+ fclose(dat);
end
end
--- /dev/null
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+#include "SLPrecangle.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);
+ plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
+
+// plhs[0] = mxCreateDoubleMatrix(1, 2, mxREAL);
+
+
+ //printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+
+ //sol = G00(p,y1,y2,x1,x2,l);
+
+ *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]));
+ *mxGetPr(plhs[1]) = FLO_perp(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+
+// printf("%f",FLO_plane(*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) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+
+
+}
+
+
+
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
//sicherheitsabfragen zu Datengroessen
if (nrhs != 7)
- mexErrMsgTxt("expected (1,x2,y1,y2,d1,d2,d3)");
+ mexErrMsgTxt("expected (x1,x2,y1,y2,d1,d2,d3)");
if (nlhs > 2)
mexErrMsgTxt("has maximal two output arguments");
mex mex_g0.cpp SLPrecangle.cpp
mex mex_G00.cpp SLPrecangle.cpp
mex mex_Fpar.cpp SLPrecangle.cpp
+mex mex_Fort.cpp SLPrecangle.cpp
mex build_A.cpp SLPrecangle.cpp
% end
-% function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
-%
-% [x1n x1w] = gauss(w,a,b);
-% [x2n x2w] = gauss(w,c,d);
-% [y1n y1w] = gauss(w,r,t);
-% [y2n y2w] = gauss(w,u,v);
-%
-% in = 0;
-% for i=1:length(y2n)
-% for j=1:length(y1n)
-% for k=1:length(x2n)
-% for l = 1:length(x1n)
-%
-% in = in + y2w(i) * y1w(j) * x2w(k) * x1w(l) * f(x1n(l),x2n(k),y1n(j),y2n(i));
-%
-% end
-% end
-% end
-% end
-%
-% end
-
function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
[x1n x1w] = gauss(w,a,b);
[x2n x2w] = gauss(w,c,d);
-
-
- d1=r-a;
- d2=u-c;
-% sum(x1w)a
- x1n=x1n-d1; x2n=x2n-d2;
-% c
-% d
-
- p = -0.5;
+ [y1n y1w] = gauss(w,r,t);
+ [y2n y2w] = gauss(w,u,v);
in = 0;
-
+ for i=1:length(y2n)
+ for j=1:length(y1n)
for k=1:length(x2n)
for l = 1:length(x1n)
- [my(1) mk(1)] = mex_G00(p,r,u,x1n(l),x2n(k),0);
- [my(2) mk(2)] = mex_G00(p,t,v,x1n(l),x2n(k),0);
- [my(3) mk(3)] = mex_G00(p,t,u,x1n(l),x2n(k),0);
- [my(4) mk(4)] = mex_G00(p,r,v,x1n(l),x2n(k),0);
-
-
-% my_sol = my(1)+my(2)-my(3)-my(4);
-
- in = in + x2w(k) * x1w(l) * (my(1)+my(2)-my(3)-my(4));
+ in = in + y2w(i) * y1w(j) * x2w(k) * x1w(l) * f(x1n(l),x2n(k),y1n(j),y2n(i));
end
end
-% in
+ end
+ end
+
+end
-end
\ No newline at end of file
+% function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
+%
+% [x1n x1w] = gauss(w,a,b);
+% [x2n x2w] = gauss(w,c,d);
+%
+%
+% d1=r-a;
+% d2=u-c;
+% % sum(x1w)a
+% x1n=x1n-d1; x2n=x2n-d2;
+% % c
+% % d
+%
+% p = -0.5;
+%
+% in = 0;
+%
+% for k=1:length(x2n)
+% for l = 1:length(x1n)
+%
+% [my(1) mk(1)] = mex_G00(p,r,u,x1n(l),x2n(k),0);
+% [my(2) mk(2)] = mex_G00(p,t,v,x1n(l),x2n(k),0);
+% [my(3) mk(3)] = mex_G00(p,t,u,x1n(l),x2n(k),0);
+% [my(4) mk(4)] = mex_G00(p,r,v,x1n(l),x2n(k),0);
+%
+%
+% % my_sol = my(1)+my(2)-my(3)-my(4);
+%
+% in = in + x2w(k) * x1w(l) * (my(1)+my(2)-my(3)-my(4));
+%
+% end
+% end
+% % in
+%
+% end
\ No newline at end of file
--- /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
+
\r
%% Laden des Gitters\r
-load exmpl_2DLShape %Energienorm sollte gegen 8.28466 konvergieren\r
+load exmpl_3DFichCube %Energienorm sollte gegen 8.28466 konvergieren\r
\r
%% Rotationen einzelnen von Elementen (zum Testen)\r
% coordinates = coordinates(:,[ 1 3 2]);\r
% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
% [coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));\r
\r
-for ref = 1:4\r
+for ref = 1:3\r
\r
%% Testen mit Quadratur... und mex_G00\r
% A = build_A2(coordinates,elements);\r
projects / bacc.git / commitdiff