else if (p == -0.5)\r
sol = asinh((y - x) / fabs(l));\r
else if (p == -1)\r
- sol = atan((y - x) / fabs(l))/2;\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
--- /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
+/* 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
+/* 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
projects / bacc.git / commitdiff