+++ /dev/null
-#include <iostream>\r
-#include <cmath>\r
-#include <cassert>\r
-#include <stdlib.h>\r
-\r
-#include "SLPrecangle.hpp"\r
-\r
-#define EPS 0.00001\r
-\r
-using namespace std;\r
-\r
-int sign(double);\r
-//double arsinh(double);\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
- * 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
- double pt = p + 1.5;\r
- double sol = 0;\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
- } 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
- cout << "warning in G00: no case for p="<< p <<" defined. Possible: [-1.5,-0.5,0.5,...]\n";\r
- //mexErrMsgTxt("no case for p defined\n");\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
- 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,\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
-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
- cout << "p must be either 0, -1/2, -1 or -3/2. (" << p << ")\n";\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
-\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
+++ /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);
-
-
-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(
- 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
+++ /dev/null
-#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 "SLPrecangle.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
- = 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
- //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
- = 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
-\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
- = 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 * 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
#include "mex.h"
#include <stdlib.h>
-#include "SLPrecangle.hpp"
+#include "slpRectangle.hpp"
#define my_PI 3.141592653589793
#define EPS 0.02
#include "mex.h"
#include <stdlib.h>
-#include "SLPrecangle.hpp"
+#include "slpRectangle.hpp"
#define my_PI 3.141592653589793
#define EPS 0.02
#include "mex.h"
#include <stdlib.h>
-#include "SLPrecangle.hpp"
+#include "slpRectangle.hpp"
#define my_PI 3.141592653589793
#define EPS 0.02
--- /dev/null
+#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
+ = apply0Int(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
+ //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
+ = apply0Int(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
+\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
+ = apply0Int(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
+ = apply0Int(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
+ = apply0Int(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
+ = apply0Int(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 * 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
#include "mex.h"
#include <stdlib.h>
-#include "SLPrecangle.hpp"
+#include "slpRectangle.hpp"
#define my_PI 3.141592653589793
#define EPS 0.02
-mex mex_g0.cpp SLPrecangle.cpp
-mex mex_G00.cpp SLPrecangle.cpp
-mex mex_Fpar.cpp SLPrecangle.cpp
-mex mex_Fort.cpp SLPrecangle.cpp
+mex mex_g0.cpp slpRectangle.cpp
+mex mex_G00.cpp slpRectangle.cpp
+mex mex_Fpar.cpp slpRectangle.cpp
+mex mex_Fort.cpp slpRectangle.cpp
-mex build_A.cpp SLPrecangle.cpp
+mex mex_build_A.cpp slpRectangle.cpp
% Diese Funktion Berechnet die Orthogonalen mit Laenge 1 über alle Flächen
% FLAG:
% w -> Laenge entspricht Flaecheninhalt
+%
% P.Schaefer
%% Parameterueberpruefung
function [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type)
+%
+% [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type)
+%
+% Verfeinert die markierten Elemente mit dem entsprechenden TYP und gibt
+% auch die F2S Beziehungen zurück. type muss von der Länge der Anzahl der
+% Elemente entsprechen und die Einträge können 1,2,3,4 sein.
+%
+% der Typ zu jedem Element entspricht dabei:
+% 1- keine Verfeinerung
+% 2- 4 neue Elemente
+% 3- 2 neue Elemente, übereinander
+% 4- 2 neue Elemente, nebeneinander
+%
+% P. Schaefer
if([1 1] == size(type))
type = repmat(type, size(elements,1),1);
--- /dev/null
+#include <iostream>\r
+#include <cmath>\r
+#include <cassert>\r
+#include <stdlib.h>\r
+\r
+#include "slpRectangle.hpp"\r
+\r
+#define EPS 0.00001\r
+\r
+using namespace std;\r
+\r
+int sign(double);\r
+//double arsinh(double);\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
+ * 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
+ double pt = p + 1.5;\r
+ double sol = 0;\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
+ } 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
+ cout << "warning in G00: no case for p="<< p <<" defined. Possible: [-1.5,-0.5,0.5,...]\n";\r
+ //mexErrMsgTxt("no case for p defined\n");\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 applyInt(\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
+ 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 apply0Int(\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
+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
+ cout << "p must be either 0, -1/2, -1 or -3/2. (" << p << ")\n";\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
+\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
--- /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);
+
+
+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 applyInt(
+ double(*)(double, double, double, double, double, double, double),
+ double, double, double, double, double, double, double, double, double,
+ double, double);
+
+// sol = quad0Int((F_par/F_ort),b,d,t,v,d1,d2,d3);
+double apply0Int(
+ double(*)(double, double, double, double, double, double, double),
+ double, double, double, double, double, double, double);
+
+
+#endif
-% function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
-% %Vierfachintegral ueber f(x1,x2,y1,y2) auf [a b]x[c d]x[r t]x[u v] mit
-% %der Simpsonregel und w Auswertungspunkten
-% s = (b-a)/w;
-% in = 0;
-% for i = a:s:b-s
-% in = in + fx2(f,i,c,d,r,t,u,v,w) + ...
-% 4*fx2(f,(2*i+s)/2,c,d,r,t,u,v,w) + ...
-% fx2(f,i+s,c,d,r,t,u,v,w);
-% end
-% in = in *s/6;
-% end
-%
-% function in = fx2(f,x1,c,d,r,t,u,v,w)
-% s = (d-c)/w;
-% in = 0;
-% for i = c:s:d-s
-% in = in + fy1(f,x1,i,r,t,u,v,w) + ...
-% 4*fy1(f,x1,(2*i+s)/2,r,t,u,v,w) + ...
-% fy1(f,x1,i+s,r,t,u,v,w);
-% end
-% in = in *s/6;
-% end
-%
-% function in = fy1(f,x1,x2,r,t,u,v,w)
-% s = (t-r)/w;
-% in = 0;
-% for i = r:s:t-s
-% in = in + fy2(f,x1,x2,i,u,v,w) + ...
-% 4*fy2(f,x1,x2,(2*i+s)/2,u,v,w) + ...
-% fy2(f,x1,x2,i+s,u,v,w);
-% end
-% in = in *s/6;
-% end
-%
-% function in = fy2(f,x1,x2,y1,u,v,w)
-% s = (v-u)/w;
-% in = 0;
-% for i = u:s:v-s
-% in = in + f(x1,x2,y1,i) + ...
-% 4*f(x1,x2,y1,(2*i+s)/2) + ...
-% f(x1,x2,y1,i+s);
-% end
-% in = in *s/6;
-% end
-
-
function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
+%
+% in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
+%
+% Berechnet das Vierfachintegral über die Funktion f(x1,x2,y1,y2) mit den
+% Grenzen [a b]x[c d]x[r t]x[u v]. Dazu wird eine Gaussquadratur über w
+% Punkte verwendet.
+%
+% P. Schaefer
[x1n x1w] = gauss(w,a,b);
[x2n x2w] = gauss(w,c,d);
for k=1:length(x2n)
for l = 1:length(x1n)
- in = in + y2w(i) * y1w(j) * x2w(k) * x1w(l) * f(x1n(l),x2n(k),y1n(j),y2n(i));
+ in = in + y2w(i) * y1w(j) * x2w(k) * x1w(l) *...
+ f(x1n(l),x2n(k),y1n(j),y2n(i));
end
end
end
end
-end
-
-% function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
-%
-% [x1n x1w] = gauss(w,a,b);
-% [x2n x2w] = gauss(w,c,d);
-%
-%
-% d1=r-a;
-% d2=u-c;
-% % sum(x1w)a
-% x1n=x1n-d1; x2n=x2n-d2;
-% % c
-% % d
-%
-% p = -0.5;
-%
-% in = 0;
-%
-% for k=1:length(x2n)
-% for l = 1:length(x1n)
-%
-% [my(1) mk(1)] = mex_G00(p,r,u,x1n(l),x2n(k),0);
-% [my(2) mk(2)] = mex_G00(p,t,v,x1n(l),x2n(k),0);
-% [my(3) mk(3)] = mex_G00(p,t,u,x1n(l),x2n(k),0);
-% [my(4) mk(4)] = mex_G00(p,r,v,x1n(l),x2n(k),0);
-%
-%
-% % my_sol = my(1)+my(2)-my(3)-my(4);
-%
-% in = in + x2w(k) * x1w(l) * (my(1)+my(2)-my(3)-my(4));
-%
-% end
-% end
-% % in
-%
-% end
\ No newline at end of file
+end
\ No newline at end of file
-% function in = surfQuad(f,a,b,c,d,v)\r
-% %Doppelintegral ueber f(x,y) auf [a b]x[c d] mit der Simpsonregel und v\r
-% %Auswertungspunkten\r
-% s = (b-a)/v;\r
-% in = 0;\r
-% for i = a:s:b-s\r
-% in = in + fy(f,i,c,d,v) + 4*fy(f,(2*i+s)/2,c,d,v) + fy(f,i+s,c,d,v);\r
-% end\r
-% in = in *s/6;\r
-% end\r
-% \r
-% function in = fy(f,x,c,d,v)\r
-% s = (d-c)/v;\r
-% in = 0;\r
-% for i = c:s:d-s\r
-% in = in + f(x,i) + 4*f(x,(2*i+s)/2) + f(x,i+s);\r
-% end\r
-% in = in *s/6;\r
-% end\r
-\r
-\r
function in = surfQuad(f,a,b,c,d,v)\r
+%\r
+% in = surfQuad(f,a,b,c,d,v)\r
+%\r
+% Berechnet das Doppelintegral über die Funktion f(x,y) mit den Grenzen\r
+% [a b]x[c d]. Dazu wird eine Gaussquadratur über v Punkte verwendet.\r
+%\r
+% P. Schaefer\r
+\r
[xn xw] = gauss(v,a,b);\r
[yn yw] = gauss(v,c,d);\r
\r
N = size(elements,1);
tic
-A = build_A(coordinates,elements);
+A = mex_build_A(coordinates,elements);
t = toc;
b = sqrt(sum(quadNorm(coordinates,elements,'w').^2,2));
ref = ref+1;\r
\r
[coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
- A_fine = build_A(coordinates_fine,elements_fine);\r
+ A_fine = mex_build_A(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
ref = ref+1;\r
\r
[coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));\r
- A_fine = build_A(coordinates_fine,elements_fine);\r
+ A_fine = mex_build_A(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
t(1) = toc;
disp ' '
tic
-A2 = build_A(coordinates,elements);
+A2 = mex_build_A(coordinates,elements);
t(2) = toc;
b = sqrt(sum(quadNorm(coordinates,elements,'w').^2,2));
projects / bacc.git / commitdiff