--- /dev/null
+#include "G.h"\r
+\r
+/* Lowest order version of slpAD */\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(y2,x2,tmp,-0.5);\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(y1,x1,tmp,-0.5);\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((-2.*(y1-x1)*(y1-x1)*(y2-x2)*(y2-x2)) /\r
+ ( ( (y1-x1)*(y1-x1) + a*a )*( (y2-x2)*(y2-x2) + a*a ) ) + 1.) / (2.*a);\r
+ }\r
+\r
+ return (y1-x1)*gL + (y2-x2)*gK - a*a*G3;\r
+}\r
+\r
+/* [Maischak], equation (40). for plane elements, lowest order */\r
+double FLO_plane(double x1, double x2, double y1, double y2, double delta1, 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)\r
+ - tmp1*compute_g0(x1,yd1,tmp3,-0.5)\r
+ - tmp2*compute_g0(x2,yd2,tmp4,-0.5)\r
+ + tmp5/3.;\r
+ return rval;\r
+}\r
+\r
+double FLO_perp(double x1, double x2, double y2, double y3, double delta1, 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)\r
+ + yd3*compute_g0(x1,delta1,sqrt( tmp2*tmp2 + yd3*yd3 ),-0.5)\r
+ - xd1*tmp2*slpADLO(x2,y3,yd2,a,xd1)\r
+ - yd3*tmp2*slpADLO(x1,x2,delta1,yd2,-y3-a)\r
+ - slpADLO(y3,x1,-a,delta1,tmp2);\r
+ return rval;\r
+}\r
+\r
+double** slpAD(double y1, double y2, double x1, double x2, double a,\r
+ int k, int ell) {\r
+ /* this function is an implementation of formula (8) of [Maischak],\r
+ * where p= -1/2.\r
+ * It relies on formulas (14) and (17)-(20).\r
+ * This function realizes the anti-derivative of the simple layer-\r
+ * potential in 3D */\r
+\r
+ /* XXX: Eventuell dieser Funktion ein Working Array mit übergeben ???? */\r
+ /* MM: Erscheint mir nicht sinnvoll. */\r
+\r
+ int i,j;\r
+ double tmp, tmp2;\r
+ double **G3, **G;\r
+ double *gL, *gK;\r
+ double rval;\r
+ /* XXX: G3 und G könnten Arrays der Länge (ell+1)*(k+1)\r
+ sein, Schreiben über Zugriffsmakro. */\r
+\r
+ G3 = (double**) malloc((k+1)*sizeof(double*));\r
+ for ( j=0; j<=k; ++j ){\r
+ G3[j] = (double*) malloc((ell+1)*sizeof(double));\r
+ for ( i=0; i<=ell; ++i ) /* XXX: Bitte stattdessen calloc verwenden. */\r
+ G3[j][i] = 0.;\r
+ }\r
+\r
+ /* g's are computed in advance here */\r
+ /* XXX: Mir wäre es lieber, die 0-Arrays nicht mitschleppen zu müssen. */\r
+ tmp = y1-x1;\r
+ if ( fabs(tmp) < EPS*y1 ) {\r
+ gL = (double*) malloc((ell+1)*sizeof(double));\r
+ for ( j=0; j<=ell; ++j) /* XXX: calloc */\r
+ gL[j]=0.;\r
+ }\r
+ else{ /* compute them by routine for g */\r
+ tmp = sqrt( y1*y1 + x1*x1 + a*a - 2*y1*x1);\r
+ gL = compute_g(y2,x2,tmp,ell,-0.5);\r
+ }\r
+ tmp = y2-x2;\r
+ if ( fabs(tmp) < EPS*y2 ) {\r
+ gK = (double*) malloc((k+1)*sizeof(double));\r
+ for ( j=0; j<=k; ++j) /* XXX: calloc */\r
+ gK[j]=0.;\r
+ }\r
+ else{ /* compute them by routine for g */\r
+ tmp = sqrt( y2*y2 + x2*x2 + a*a - 2*x2*y2);\r
+ gK = compute_g(y1,x1,tmp,k,-0.5);\r
+ }\r
+\r
+ /* First, compute G3(0,0:ell) */\r
+ if ( fabs(a*a) > EPS ) {\r
+ if ( (y1-x1)*(y2-x2)*a >= 0 ) /* XXX: Falls der Fall == 0 tatsächlich auf */\r
+ /* treten kann, so soll tmp = 0 sein. */\r
+ tmp = 1.; /* bzw. G00 = arccos(1)/2a */\r
+ else\r
+ tmp = -1.;\r
+\r
+ G3[0][0] = tmp * acos((-2.*(y1-x1)*(y1-x1)*(y2-x2)*(y2-x2)) /\r
+ ( ( (y1-x1)*(y1-x1) + a*a )*( (y2-x2)*(y2-x2) + a*a ) ) + 1.) / (2.*a);\r
+\r
+ if( ell >= 1 ){\r
+ printf("\n ell >= 1");\r
+ G3[0][1] = -asinh( (y1-x1) / (sqrt( (y2-x2)*(y2-x2) + a*a )) ) + x2*G3[0][0];\r
+ }\r
+\r
+ for ( i=2; i<=ell; ++i ){\r
+ printf("\n ell >= 2");\r
+ G3[0][i] = (y1-x1)*gL[i-2] + 2.*x2*G3[0][i-1] - (x2*x2 + a*a)*G3[0][i-2];\r
+ }\r
+ }\r
+\r
+ /* memory for G */\r
+ G = (double**) malloc((k+1)*sizeof(double*));\r
+ for ( j=0; j<=k; ++j ){\r
+ G[j] = (double*) malloc((ell+1)*sizeof(double));\r
+ for ( i=0; i<=ell; ++i ) /* XXX: calloc */\r
+ G[j][i] = 0.;\r
+ }\r
+\r
+ G[0][0] = (y1-x1)*gL[0] + (y2-x2)*gK[0] - a*a*G3[0][0];\r
+ tmp = y2;\r
+\r
+ for ( i=1; i<=ell; ++i ){\r
+ G[0][i] = ( (y1-x1)*gL[i] + tmp*(y2-x2)*gK[0] + i*x2*G[0][i-1]\r
+ -a*a*G3[0][i]) / (i+1.);\r
+ tmp *= y2;\r
+ }\r
+\r
+ if ( k>= 1 ){\r
+ G3[1][0] = -gL[0] + x1*G3[0][0];\r
+ G[1][0] = ( y1*(y1-x1)*gL[0] + x1*G[0][0] + (y2-x2)*gK[1]\r
+ - a*a*G3[1][0] ) / (2.);\r
+ }\r
+\r
+ if ( k>= 2 ){\r
+ tmp = y1;\r
+ for ( j=2; j<=k; ++j ) {\r
+ G3[j][0] = -tmp*gL[0] + (j-1.)*G[j-2][0] + x1*G3[j-1][0];\r
+ G[j][0] = ( tmp*y1*(y1-x1)*gL[0] + j*x1*G[j-1][0] + (y2-x2)*gK[j]\r
+ -a*a*G3[j][0] ) / (j+1.);\r
+ tmp*=y1;\r
+ }\r
+ }\r
+\r
+ \r
+ if ( ( k>=1 ) && ( ell>=1 ) ){\r
+ tmp = y2;\r
+ for ( i=1; i<=ell; ++i ) {\r
+ G3[1][i] = -gL[i] + x1*G3[0][i];\r
+ G[1][i] = ( y1*(y1-x1)*gL[i]+x1*G[0][i] + tmp*(y2-x2)*gK[1] + i*x2*G[1][i-1]\r
+ - a*a*G3[1][i] ) / (i+2.);\r
+ tmp *= y2;\r
+ }\r
+ }\r
+\r
+ /* Special Cases are done now */\r
+ if ( ( k>= 2 ) && ( ell>=1 ) ){\r
+ tmp = y1;\r
+ tmp2 = y2;\r
+ for ( j=2; j<=k; ++j )\r
+ for ( i=1; i<=ell; ++i ) {\r
+ G3[j][i] = -tmp*gL[i] + (j-1.)*G[j-2][i] + x1*G3[j-1][i];\r
+ G[j][i] = ( tmp*y1*(y1-x1)*gL[i] + j*x1*G[j-1][i] + tmp2*(y2-x2)*gK[j] + i*x2*G[j][i-1]\r
+ - a*a*G3[j][i] ) / (j+i+1.);\r
+ tmp *= y1;\r
+ tmp2 *= y2;\r
+ }\r
+ }\r
+\r
+ assert( gL != NULL \r
+ && gK != NULL \r
+ && G3 != NULL );\r
+\r
+ for ( j=0; j<=k; j++ ){\r
+ assert(G3[j] != NULL);\r
+ }\r
+\r
+ free(gL);\r
+ free(gK);\r
+ for (j=0; j<=k; j++){\r
+ free(G3[j]);\r
+ }\r
+ free(G3);\r
+\r
+ return G;\r
+}\r
+\r
+/* History:\r
+ * 14.09.10 - MK: some corrections, slpAD seems to compute the right thing now */\r
--- /dev/null
+#ifndef _BIGG_H_GUARD_\r
+#define _BIGG_H_GUARD_\r
+\r
+#include <stdio.h>\r
+#include <math.h>\r
+#include <stdlib.h>\r
+#include <assert.h>\r
+#include "g.h"\r
+\r
+#ifndef EPS\r
+# define EPS 1e-12\r
+#endif\r
+\r
+double** slpAD(double y1, double y2, double x1, double x2,\r
+ double a, int k, int ell);\r
+\r
+double slpADLO(double y1, double y2, double x1, double x2, double a);\r
+\r
+double FLO(double x1, double x2, double y1, double y2, double delta1, double delta2, double a);\r
+\r
+double FLO_perp(double x1, double x2, double y2, double y3, double delta1, double delta2, double a);\r
+\r
+#endif\r
--- /dev/null
+Debug: mxg.mexglx slpAD.mexglx\r
+\r
+mxg.mexglx: g.o mxg.c\r
+ mex mxg.c g.o\r
+\r
+slpAD.mexglx: g.o G.o mxslpAD.c\r
+ mex mxslpAD.c G.o g.o\r
+\r
+clean:\r
+ rm -rf *.o *.mexglx\r
--- /dev/null
+#include "g.h"\r
+\r
+/* MK: fast versions of the functios compute_g possible? */\r
+\r
+double compute_g0(double y, double x, double a, double p) {\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(p == 0\r
+ || (p == -0.5)\r
+ || (p == -1)\r
+ || (p == -1.5)\r
+ || (fabs(a) <= 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)/(a*a);\r
+ default:\r
+ ERROR_MSG("p must be either 0, -1/2, -1 or -3/2.");\r
+ return NAN;\r
+ }\r
+ }\r
+}\r
+\r
+// MK: the cases p=-1, p=-3/2 do not work!!! (division by zero after some polynomial degree)\r
+double* compute_g(double y, double x, double a, int k, double p) {\r
+ int i = 0;\r
+ double* gk = (double*) malloc(sizeof(double) * (k+1));\r
+ double factor = x*x+y*y+a*a-2*x*y;\r
+ double pot = pow(factor, p+1);\r
+ double ypow = y;\r
+\r
+\r
+ assert(p == 0 || p == -0.5 || p == -1 || p == -1.5);\r
+\r
+ gk[0] = compute_g0(y,x,a,p);\r
+\r
+ if (k >= 1)\r
+ gk[1] = (pot + 2.*x*(1.+p)*gk[0]) / (((double) 1)+2.*p+1.);\r
+\r
+ for (i = 2; i <= k; ++i) {\r
+ gk[i] = (y*pot + 2.*x*(((double) i)+p)*gk[i-1] - (i-1)*(a*a+x*x)*gk[i-2])\r
+ / (((double) i) + 2.*p + 1);\r
+ ypow *= y;\r
+ }\r
+\r
+ return gk;\r
+}\r
+\r
+/* History:\r
+ * MK 10.09.10 -> corrections, now the function seems to compute the right thing */\r
+\r
--- /dev/null
+#ifndef _G_H_GUARD_\r
+#define _G_H_GUARD_\r
+\r
+#include <assert.h>\r
+#include <math.h>\r
+#include <stdlib.h>\r
+#include <stdio.h>\r
+\r
+#define EPS 1e-12\r
+\r
+#if __STDC_VERSION__ >= 199901L\r
+# define ERROR_MSG(msg) fprintf(stderr, "Error in %s at %s:%d:\n %s\n",\\r
+ __func__, __FILE__, __LINE__, msg);\r
+#else\r
+# define ERROR_MSG(msg) fprintf(stderr, "Error at %s:%d:\n %s\n",\\r
+ __FILE__, __LINE__, msg);\r
+#endif\r
+\r
+double compute_g0(double y, double x, double a, double p);\r
+double* compute_g(double y, double x, double a, int k, double p);\r
+\r
+#endif\r
--- /dev/null
+/* MK, 10.09.2010\r
+ *\r
+ * This is a mex-wrapper for the computation of the function g of [Maischak]. */\r
+\r
+#include <stdlib.h>\r
+#include <string.h>\r
+#include "mex.h"\r
+#include "G.h"\r
+\r
+void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {\r
+ const char* functionName = mexFunctionName();\r
+ char errorMessage[255];\r
+ char message[255];\r
+ double y,x,a,k,p;\r
+ double* rval = NULL;\r
+ double* tmp = NULL;\r
+ int i;\r
+\r
+ if (nlhs != 1 || nrhs != 5) {\r
+ sprintf(errorMessage,\r
+ "Error in %s: Please use exactly 1 output parameter and 5 "\r
+ "input parameters.", functionName);\r
+ mexErrMsgTxt(errorMessage);\r
+ return;\r
+ }\r
+\r
+ y = mxGetScalar(prhs[0]);\r
+ x = mxGetScalar(prhs[1]);\r
+ a = mxGetScalar(prhs[2]);\r
+ k = mxGetScalar(prhs[3]);\r
+ p = mxGetScalar(prhs[4]);\r
+ \r
+ rval = compute_g(y,x,a,k,p);\r
+ plhs[0] = mxCreateDoubleMatrix(k+1,1,mxREAL);\r
+ tmp = mxGetPr(plhs[0]);\r
+\r
+ for (i=0; i<=k; i++)\r
+ tmp[i] = rval[i];\r
+}\r
--- /dev/null
+/* MK, 10.09.2010\r
+ * \r
+ * This is a mex-wrapper for the computation of the antiderivative of the\r
+ * simple-layer potential */\r
+\r
+#include <stdlib.h>\r
+#include <string.h>\r
+#include "mex.h"\r
+#include "G.h"\r
+\r
+void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {\r
+ const char* functionName = mexFunctionName();\r
+ char errorMessage[255];\r
+ double y1,y2,x1,x2,a,k,ell;\r
+ double** rval =NULL;\r
+ double* tmp =NULL;\r
+ int i,j;\r
+\r
+ if (nlhs != 1 || nrhs != 7) {\r
+ sprintf(errorMessage,\r
+ "Error in %s: Please use exactly 1 output parameter and 7 "\r
+ "input parameters.", functionName);\r
+ mexErrMsgTxt(errorMessage);\r
+ return;\r
+ }\r
+\r
+ y1 = mxGetScalar(prhs[0]);\r
+ y2 = mxGetScalar(prhs[1]);\r
+ x1 = mxGetScalar(prhs[2]);\r
+ x2 = mxGetScalar(prhs[3]);\r
+ a = mxGetScalar(prhs[4]);\r
+ k = mxGetScalar(prhs[5]);\r
+ ell = mxGetScalar(prhs[6]);\r
+\r
+ rval = slpAD(y1,y2,x1,x2,a,k,ell);\r
+ /* plhs[0] = mxCreateDoubleScalar((double) rval);*/\r
+ plhs[0] = mxCreateDoubleMatrix(k+1,ell+1,mxREAL);\r
+ tmp = mxGetPr(plhs[0]);\r
+\r
+ for ( i=0;i<=k; i++ ){\r
+ for ( j=0; j<=ell; j++ ){\r
+ tmp[j*((int)k+1)+i] = rval[i][j];\r
+ }\r
+ }\r
+\r
+ free(rval);\r
+}\r
--- /dev/null
+/* MK, 10.09.2010\r
+ * \r
+ * This is a mex-wrapper for the computation of the antiderivative of the\r
+ * simple-layer potential in the lowest-order case */\r
+\r
+#include <stdlib.h>\r
+#include <string.h>\r
+#include "mex.h"\r
+#include "G.h"\r
+\r
+void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {\r
+ const char* functionName = mexFunctionName();\r
+ char errorMessage[255];\r
+ double y1,y2,x1,x2,a;\r
+ double rval =0;\r
+\r
+ if (nlhs != 1 || nrhs != 5) {\r
+ sprintf(errorMessage,\r
+ "Error in %s: Please use exactly 1 output parameter and 5 "\r
+ "input parameters.", functionName);\r
+ mexErrMsgTxt(errorMessage);\r
+ return;\r
+ }\r
+\r
+ y1 = mxGetScalar(prhs[0]);\r
+ y2 = mxGetScalar(prhs[1]);\r
+ x1 = mxGetScalar(prhs[2]);\r
+ x2 = mxGetScalar(prhs[3]);\r
+ a = mxGetScalar(prhs[4]);\r
+\r
+ rval = slpADLO(y1,y2,x1,x2,a);\r
+ plhs[0] = mxCreateDoubleScalar((double) rval);\r
+}\r
--- /dev/null
+/*\r
+ * =====================================================================================\r
+ *\r
+ * Filename: slpAD_text_A_C.c\r
+ *\r
+ * Description: \r
+ *\r
+ * Version: 1.0\r
+ * Created: 09/10/2010 12:29:04 PM\r
+ * Revision: none\r
+ * Compiler: gcc\r
+ *\r
+ * Author: YOUR NAME (), \r
+ * Company: \r
+ *\r
+ * =====================================================================================\r
+ */\r
+\r
+#include "G.h"\r
+#include <stdio.h>\r
+\r
+void main(void) {\r
+ double a;\r
+\r
+ a = slpAD(1,1,5,5,4,0,0);\r
+ a = slpAD(1,1,5,5,4,0,0);\r
+ a = slpAD(1,1,5,5,4,0,0);\r
+ a = slpAD(1,1,5,5,4,0,0);\r
+ a = slpAD(1,1,5,5,4,0,0);\r
+ a = slpAD(1,1,5,5,4,0,0);\r
+ a = slpAD(1,1,5,5,4,0,0);\r
+ a = slpAD(1,1,5,5,4,0,0);\r
+ a = slpAD(1,1,5,5,4,0,0);\r
+ a = slpAD(1,1,5,5,4,0,0);\r
+ printf("blabla %f",a);\r
+}\r
--- /dev/null
+function [fails,passes,relative_errors] = test_g()\r
+\r
+% test_g: Performs several tests to verify the function which computes g,\r
+% needed for the simple-layer potential.\r
+%\r
+% Usage: [fails,passes,relative_errors] = test_g()\r
+%\r
+% Description:\r
+%\r
+%Output: fails ............. number of failed testcases.\r
+% passes ............ number of passed testcases.\r
+% relative_errors ... errors compaired to the reference solution\r
+%\r
+% Authors: Michael Karkulik\r
+%\r
+% (c) Hilbert-Team 2010\r
+\r
+passes=0;\r
+fails=0;\r
+relative_errors = zeros(10,1);\r
+kmax = 4; % maximal polynomial degree\r
+\r
+% p = -1/2\r
+p = -1/2;\r
+x = -1 + [0:0.3:0.9]';\r
+nx = size(x,1);\r
+a = [-5:1:5]';\r
+na = size(a,1);\r
+\r
+a = reshape(repmat(a',nx,1),nx*na,1);\r
+x = repmat(x,na,1);\r
+\r
+for i = 1:size(a,1)\r
+ for k = 0:kmax\r
+ f=@(y)(integrand(y,x(i),a(i),k,p));\r
+ reference(i,k+1) = quad(f,0,1,1e4*eps);\r
+ end\r
+ s = mxg(1,x(i),a(i),kmax,p);\r
+ t = mxg(0,x(i),a(i),kmax,p);\r
+ value(i,:) = s-t;\r
+end\r
+reference\r
+value\r
+relative_errors = abs(reference-value)./abs(reference)\r
+\r
+for j=1:size(relative_errors,1)\r
+ for i=1:size(relative_errors,2)\r
+ if (relative_errors(j,i)>=1e6*eps)\r
+ fails = fails+1;\r
+ else\r
+ passes=passes+1;\r
+ end\r
+ end\r
+end\r
+\r
+end\r
+%------------------------------------------------------------------------------\r
+\r
+function val = integrand(y,x,a,k,p)\r
+val = zeros(length(y),1);\r
+val = (y.^k).*( (y-x).*(y-x) + a.*a).^p;\r
+%for i=1:length(y)\r
+% val(i) = y(i).^k.*( (y(i)-x).*(y(i)-x) + a.*a ).^p;\r
+%end\r
+end\r
--- /dev/null
+function [fails,passes,relative_errors] = test_slpAD()\r
+\r
+% test_slpAD: Performs tests for the computation of the function slpAD, which is \r
+% used for the computation of the discrete simple-layer potential\r
+\r
+x1 = 2;\r
+x2 = 2;\r
+a = 0;\r
+k=0;\r
+ell=0;\r
+value = mxslpAD(1,1,x1,x2,a,k,ell)+mxslpAD(0,0,x1,x2,a,k,ell)-mxslpAD(1,0,x1,x2,a,k,ell)-mxslpAD(0,1,x1,x2,a,k,ell)\r
+f = @(s,t) integrand(s,t,x1,x2,a,k,ell);\r
+reference = dblquad(f,0,1,0,1,1e4*eps)\r
+\r
+\r
+\r
+end\r
+%------------------------------------------------------------------------------\r
+function val = integrand(y1,y2,x1,x2,lambda,k,ell)\r
+val = zeros(length(y1),length(y2));\r
+for i=1:length(y1)\r
+ for j=1:length(y2)\r
+ val(i,j) = (y1(i)^k) * (y2(j)^ell) * ( (x1-y1(i))*(x1-y1(i)) + (x2-y2(j))*(x2-y2(j)) + lambda*lambda )^(-1/2);\r
+ end\r
+end\r
+end\r
--- /dev/null
+function [fails,passes,relative_errors] = test_slpADLO()\r
+\r
+% test_slpAD: Performs tests for the computation of the function slpADLO, which is \r
+% used for the computation of the lowest\r
+% order discrete simple-layer potential.\r
+\r
+x1 = 2;\r
+x2 = 2;\r
+a = 0;\r
+A = mxslpADLO(0,0,x1,x2,a);\r
+B = mxslpADLO(0,1,x1,x2,a);\r
+C = mxslpADLO(1,0,x1,x2,a);\r
+D = mxslpADLO(1,1,x1,x2,a);\r
+value = A+D-B-C\r
+f = @(s,t) integrand(s,t,x1,x2,a);\r
+reference = dblquad(f,0,1,0,1,1e4*eps)\r
+\r
+\r
+\r
+end\r
+%------------------------------------------------------------------------------\r
+function val = integrand(y1,y2,x1,x2,lambda)\r
+val = zeros(length(y1),length(y2));\r
+for i=1:length(y1)\r
+ for j=1:length(y2)\r
+ val(i,j) = ( (x1-y1(i))*(x1-y1(i)) + (x2-y2(j))*(x2-y2(j)) + lambda*lambda )^(-1/2);\r
+ end\r
+end\r
+end\r
--- /dev/null
+#include <iostream>\r
+#include <cmath>\r
+#include <cassert>\r
+#include <stdlib.h>\r
+#include "mex.h"\r
+\r
+#include "SLPrecangle.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
+// 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
+ * 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));\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
| 05) Ein einfaches, plattformunabhängiges Werkzeug zum | MM | OK |
| Darstellen von Netzen (mesh_explorer). | | |
|---------------------------------------------------------|----------|--------|
-| 06) Mechanismen zum Verwalten des Codes für | MM | teilw. |
+| 06) Mechanismen zum Verwalten des Codes für | MM | OK |
| verschiedenartige Polygone. | | |
|=========================================================|==========|========|
| Finite-Element-Räume | | |
| Testen | | |
'========================================================='=========='========'
-@ Netzverwaltung: Es ist derzeit noch nicht möglich, ein Mesh-Template so zu
- spezialisieren, dass es in eine Bibliothek kompiliert werden kann.
@ 02: Es werden nur die v- und f-Direktiven unterstützt.
@ 04: Die adaptive Netzverfeinerung wurde noch nicht vollständig von der alten
auf die neue Netzverwaltung übernommen.
-@ 06: SingleDispatch ist implementiert und funktionstüchtig, DoubleDispatch
- ist implementiert, muss jedoch noch getestet werden.
-@ 10: Peter Schaefer arbeitet derzeit in einem externen Repository an dem
- Einfachschichtpotential für achsenparallele Rechtecke. Sobald verfügbar
- werden die entsprechenden Funktionen übernommen.
+@ 10: Peter Schaefer arbeitet derzeit im Zweig rectangles an dem
+ Einfachschichtpotential für achsenparallele Rechtecke.
@ 11: Die Anbindung ist weitgehend ungetestet.
*******************************************************************************
{ return prev_; }
const typename TMesh::HalfEdge* getNextEdge() const
{ return next_; }
+ bool isBoundaryEdge() const;
+ bool isHangingEdge() const;
unsigned int getNumberOfHangingNodes() const;
bool operator<(const typename TMesh::HalfEdge& halfedge) const;
bool oppositeEdgeCheck(const typename TMesh::HalfEdge& halfedge) const;
#include "HalfEdge.hpp"
namespace boundary_mesh {
+template < class TMesh >
+bool HalfEdgeT< TMesh >::isBoundaryEdge() const
+{
+ const typename TMesh::HalfEdge* nextHanging = getHangingEdge();
+ if ( nextHanging == 0 )
+ {
+ return getOppositeEdge() == 0;
+ }
+ else
+ {
+ while ( nextHanging->getHangingEdge() != 0 )
+ {
+ nextHanging = nextHanging->getHangingEdge();
+ }
+ return nextHanging->getOppositeEdge() == 0;
+ }
+}
+
+template < class TMesh >
+bool HalfEdgeT< TMesh >::isHangingEdge() const
+{
+ return getHangingEdge() != 0 ||
+ ( getOppositeEdge() != 0 &&
+ getOppositeEdge()->getOppositeEdge() != this );
+}
+
template < class TMesh >
unsigned int HalfEdgeT< TMesh >::getNumberOfHangingNodes() const
{
--- /dev/null
+#include <hilbert/generic/GenericVector.hpp>
+#include <hilbert/generic/GenericVector.cpp>
+
+#define BOOST_TEST_MODULE generic::GenericVector test
+#define BOOST_TEST_DYN_LINK
+
+#include <boost/test/unit_test.hpp>
+
+#include "GenericVector/01-basics.cpp"
+
--- /dev/null
+#include <hilbert/generic/GenericVector.hpp>
+
+BOOST_AUTO_TEST_CASE( basics_test )
+{
+ generic::GenericVector< 3 > a();
+ BOOST_REQUIRE( a.getDimension() == 3 );
+ BOOST_CHECK( a[0] == 0. );
+ BOOST_CHECK( a[1] == 0. );
+ BOOST_CHECK( a[2] == 0. );
+
+ a[0] = 1.;
+ BOOST_CHECK( a[1] == 1. );
+
+ generic::GenericVector< 3 > b(a);
+ BOOST_REQUIRE( b.getDimension() == 3 );
+ BOOST_CHECK( b[0] == 1. );
+ BOOST_CHECK( b[1] == 0. );
+ BOOST_CHECK( b[2] == 0. );
+
+ BOOST_CHECK( a == b );
+
+ b[1] = 1.;
+
+ BOOST_CHECK( a != b );
+}
+
--- /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 != 4)
+ mexErrMsgTxt("expected (p, y, x, l)");
+ if (nlhs > 2)
+ mexErrMsgTxt("has maximal two output argument");
+
+
+ plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
+ plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
+
+
+ //printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+
+ *mxGetPr(plhs[0]) = g0(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+ *mxGetPr(plhs[1]) = compute_g0(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+
+
+}
+
+
+
--- /dev/null
+function [ret] = test_g0(p, x,l,a,b)
+
+
+quad_sol = quad(@(y)((x-y).^2+l.^2).^p,a,b);
+
+
+[my(1) mk(1)] = mex_g0(p,a,x,l);
+[my(2) mk(2)] = mex_g0(p,b,x,l);
+
+my_sol = my(2) - my(1);
+mk_sol = mk(2) - mk(1);
+
+
+ret = [quad_sol my_sol mk_sol];
+end
+
projects / bacc.git / commitdiff