+++ /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