From 91967c6dc0803c8308eaa2e6f1a57a839ca3b51c Mon Sep 17 00:00:00 2001 From: treecity Date: Wed, 11 May 2011 18:30:04 +0000 Subject: [PATCH] =?utf8?q?[src]=20MK=20hinzugef=C3=BCgt=20[src]=20g0=20in?= =?utf8?q?=20Fall=20p=3D-1=20/lambda=20statt=20/2?= MIME-Version: 1.0 Content-Type: text/plain; charset=utf8 Content-Transfer-Encoding: 8bit git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@24 26120e32-c555-405d-b3e1-1f783fb42516 --- src/SLPrecangle.cpp | 2 +- src/maik/Makefile | 10 ++ src/maik/g.c | 215 ++++++++++++++++++++++++++++++++++++++ src/maik/g.h | 23 ++++ src/maik/mxg.c | 39 +++++++ src/maik/mxslpAD.c | 47 +++++++++ src/maik/mxslpADLO.c | 33 ++++++ src/maik/slpAD_test_A_C.c | 36 +++++++ src/maik/test_g.m | 65 ++++++++++++ src/maik/test_slpAD.m | 26 +++++ src/maik/test_slpADLO.m | 29 +++++ 11 files changed, 524 insertions(+), 1 deletion(-) create mode 100644 src/maik/Makefile create mode 100644 src/maik/g.c create mode 100644 src/maik/g.h create mode 100644 src/maik/mxg.c create mode 100644 src/maik/mxslpAD.c create mode 100644 src/maik/mxslpADLO.c create mode 100644 src/maik/slpAD_test_A_C.c create mode 100644 src/maik/test_g.m create mode 100644 src/maik/test_slpAD.m create mode 100644 src/maik/test_slpADLO.m diff --git a/src/SLPrecangle.cpp b/src/SLPrecangle.cpp index d5c9e11..4a67229 100644 --- a/src/SLPrecangle.cpp +++ b/src/SLPrecangle.cpp @@ -182,7 +182,7 @@ double g0(double p, double y, double x, double l) { else if (p == -0.5) sol = asinh((y - x) / fabs(l)); else if (p == -1) - sol = atan((y - x) / fabs(l))/2; + sol = atan((y - x) / fabs(l))/l; else if (p == -1.5) sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l)); else diff --git a/src/maik/Makefile b/src/maik/Makefile new file mode 100644 index 0000000..cee0636 --- /dev/null +++ b/src/maik/Makefile @@ -0,0 +1,10 @@ +Debug: mxg.mexglx slpAD.mexglx + +mxg.mexglx: g.o mxg.c + mex mxg.c g.o + +slpAD.mexglx: g.o G.o mxslpAD.c + mex mxslpAD.c G.o g.o + +clean: + rm -rf *.o *.mexglx diff --git a/src/maik/g.c b/src/maik/g.c new file mode 100644 index 0000000..e0b6699 --- /dev/null +++ b/src/maik/g.c @@ -0,0 +1,215 @@ +#include "G.h" + +/* Lowest order version of slpAD */ +double slpADLO(double y1, double y2, double x1, double x2, double a){ + double G3 = 0; + double gL = 0; + double gK = 0; + double tmp; + + tmp = y1-x1; + if ( fabs(tmp) >= EPS*y1 ) { + tmp = sqrt( y1*y1 + x1*x1 + a*a - 2*y1*x1); + gL = compute_g0(y2,x2,tmp,-0.5); + } + tmp = y2-x2; + if ( fabs(tmp) >= EPS*y2 ) { + tmp = sqrt( y2*y2 + x2*x2 + a*a - 2*y2*x2); + gK = compute_g0(y1,x1,tmp,-0.5); + } + if ( fabs(a*a) > EPS ) { + if ( (y1-x1)*(y2-x2)*a >= 0 ) + tmp = 1.; + else + tmp = -1.; + + G3 = tmp * acos((-2.*(y1-x1)*(y1-x1)*(y2-x2)*(y2-x2)) / + ( ( (y1-x1)*(y1-x1) + a*a )*( (y2-x2)*(y2-x2) + a*a ) ) + 1.) / (2.*a); + } + + return (y1-x1)*gL + (y2-x2)*gK - a*a*G3; +} + +/* [Maischak], equation (40). for plane elements, lowest order */ +double FLO_plane(double x1, double x2, double y1, double y2, double delta1, double delta2, double a){ + double yd1 = y1+delta1; + double yd2 = y2+delta2; + double tmp1 = x1-y1-delta1; + double tmp2 = x2-y2-delta2; + double tmp3 = sqrt( tmp2*tmp2 + a*a ); + double tmp4 = sqrt( tmp1*tmp1 + a*a ); + double tmp5 = pow(tmp1*tmp1 + tmp2*tmp2 + a*a,3./2.); + double rval = 0; + + rval = tmp1*tmp2*slpADLO(x1,x2,yd1,yd2,a) + - tmp1*compute_g0(x1,yd1,tmp3,-0.5) + - tmp2*compute_g0(x2,yd2,tmp4,-0.5) + + tmp5/3.; + return rval; +} + +double FLO_perp(double x1, double x2, double y2, double y3, double delta1, double delta2, double a){ + double xd1 = x1-delta1; + double xd2 = x2-delta2; + double yd2 = y2-delta2; + double yd3 = y3-a; + double tmp2 = x2-y2-delta2; + double rval = 0; + + rval = xd1*compute_g0(y3,a,sqrt( xd1*xd1 + tmp2*tmp2 ),-0.5) + + yd3*compute_g0(x1,delta1,sqrt( tmp2*tmp2 + yd3*yd3 ),-0.5) + - xd1*tmp2*slpADLO(x2,y3,yd2,a,xd1) + - yd3*tmp2*slpADLO(x1,x2,delta1,yd2,-y3-a) + - slpADLO(y3,x1,-a,delta1,tmp2); + return rval; +} + +double** slpAD(double y1, double y2, double x1, double x2, double a, + int k, int ell) { + /* this function is an implementation of formula (8) of [Maischak], + * where p= -1/2. + * It relies on formulas (14) and (17)-(20). + * This function realizes the anti-derivative of the simple layer- + * potential in 3D */ + + /* XXX: Eventuell dieser Funktion ein Working Array mit übergeben ???? */ + /* MM: Erscheint mir nicht sinnvoll. */ + + int i,j; + double tmp, tmp2; + double **G3, **G; + double *gL, *gK; + double rval; + /* XXX: G3 und G könnten Arrays der Länge (ell+1)*(k+1) + sein, Schreiben über Zugriffsmakro. */ + + G3 = (double**) malloc((k+1)*sizeof(double*)); + for ( j=0; j<=k; ++j ){ + G3[j] = (double*) malloc((ell+1)*sizeof(double)); + for ( i=0; i<=ell; ++i ) /* XXX: Bitte stattdessen calloc verwenden. */ + G3[j][i] = 0.; + } + + /* g's are computed in advance here */ + /* XXX: Mir wäre es lieber, die 0-Arrays nicht mitschleppen zu müssen. */ + tmp = y1-x1; + if ( fabs(tmp) < EPS*y1 ) { + gL = (double*) malloc((ell+1)*sizeof(double)); + for ( j=0; j<=ell; ++j) /* XXX: calloc */ + gL[j]=0.; + } + else{ /* compute them by routine for g */ + tmp = sqrt( y1*y1 + x1*x1 + a*a - 2*y1*x1); + gL = compute_g(y2,x2,tmp,ell,-0.5); + } + tmp = y2-x2; + if ( fabs(tmp) < EPS*y2 ) { + gK = (double*) malloc((k+1)*sizeof(double)); + for ( j=0; j<=k; ++j) /* XXX: calloc */ + gK[j]=0.; + } + else{ /* compute them by routine for g */ + tmp = sqrt( y2*y2 + x2*x2 + a*a - 2*x2*y2); + gK = compute_g(y1,x1,tmp,k,-0.5); + } + + /* First, compute G3(0,0:ell) */ + if ( fabs(a*a) > EPS ) { + if ( (y1-x1)*(y2-x2)*a >= 0 ) /* XXX: Falls der Fall == 0 tatsächlich auf */ + /* treten kann, so soll tmp = 0 sein. */ + tmp = 1.; /* bzw. G00 = arccos(1)/2a */ + else + tmp = -1.; + + G3[0][0] = tmp * acos((-2.*(y1-x1)*(y1-x1)*(y2-x2)*(y2-x2)) / + ( ( (y1-x1)*(y1-x1) + a*a )*( (y2-x2)*(y2-x2) + a*a ) ) + 1.) / (2.*a); + + if( ell >= 1 ){ + printf("\n ell >= 1"); + G3[0][1] = -asinh( (y1-x1) / (sqrt( (y2-x2)*(y2-x2) + a*a )) ) + x2*G3[0][0]; + } + + for ( i=2; i<=ell; ++i ){ + printf("\n ell >= 2"); + G3[0][i] = (y1-x1)*gL[i-2] + 2.*x2*G3[0][i-1] - (x2*x2 + a*a)*G3[0][i-2]; + } + } + + /* memory for G */ + G = (double**) malloc((k+1)*sizeof(double*)); + for ( j=0; j<=k; ++j ){ + G[j] = (double*) malloc((ell+1)*sizeof(double)); + for ( i=0; i<=ell; ++i ) /* XXX: calloc */ + G[j][i] = 0.; + } + + G[0][0] = (y1-x1)*gL[0] + (y2-x2)*gK[0] - a*a*G3[0][0]; + tmp = y2; + + for ( i=1; i<=ell; ++i ){ + G[0][i] = ( (y1-x1)*gL[i] + tmp*(y2-x2)*gK[0] + i*x2*G[0][i-1] + -a*a*G3[0][i]) / (i+1.); + tmp *= y2; + } + + if ( k>= 1 ){ + G3[1][0] = -gL[0] + x1*G3[0][0]; + G[1][0] = ( y1*(y1-x1)*gL[0] + x1*G[0][0] + (y2-x2)*gK[1] + - a*a*G3[1][0] ) / (2.); + } + + if ( k>= 2 ){ + tmp = y1; + for ( j=2; j<=k; ++j ) { + G3[j][0] = -tmp*gL[0] + (j-1.)*G[j-2][0] + x1*G3[j-1][0]; + G[j][0] = ( tmp*y1*(y1-x1)*gL[0] + j*x1*G[j-1][0] + (y2-x2)*gK[j] + -a*a*G3[j][0] ) / (j+1.); + tmp*=y1; + } + } + + + if ( ( k>=1 ) && ( ell>=1 ) ){ + tmp = y2; + for ( i=1; i<=ell; ++i ) { + G3[1][i] = -gL[i] + x1*G3[0][i]; + G[1][i] = ( y1*(y1-x1)*gL[i]+x1*G[0][i] + tmp*(y2-x2)*gK[1] + i*x2*G[1][i-1] + - a*a*G3[1][i] ) / (i+2.); + tmp *= y2; + } + } + + /* Special Cases are done now */ + if ( ( k>= 2 ) && ( ell>=1 ) ){ + tmp = y1; + tmp2 = y2; + for ( j=2; j<=k; ++j ) + for ( i=1; i<=ell; ++i ) { + G3[j][i] = -tmp*gL[i] + (j-1.)*G[j-2][i] + x1*G3[j-1][i]; + 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] + - a*a*G3[j][i] ) / (j+i+1.); + tmp *= y1; + tmp2 *= y2; + } + } + + assert( gL != NULL + && gK != NULL + && G3 != NULL ); + + for ( j=0; j<=k; j++ ){ + assert(G3[j] != NULL); + } + + free(gL); + free(gK); + for (j=0; j<=k; j++){ + free(G3[j]); + } + free(G3); + + return G; +} + +/* History: + * 14.09.10 - MK: some corrections, slpAD seems to compute the right thing now */ diff --git a/src/maik/g.h b/src/maik/g.h new file mode 100644 index 0000000..aec9e09 --- /dev/null +++ b/src/maik/g.h @@ -0,0 +1,23 @@ +#ifndef _BIGG_H_GUARD_ +#define _BIGG_H_GUARD_ + +#include +#include +#include +#include +#include "g.h" + +#ifndef EPS +# define EPS 1e-12 +#endif + +double** slpAD(double y1, double y2, double x1, double x2, + double a, int k, int ell); + +double slpADLO(double y1, double y2, double x1, double x2, double a); + +double FLO(double x1, double x2, double y1, double y2, double delta1, double delta2, double a); + +double FLO_perp(double x1, double x2, double y2, double y3, double delta1, double delta2, double a); + +#endif diff --git a/src/maik/mxg.c b/src/maik/mxg.c new file mode 100644 index 0000000..c104919 --- /dev/null +++ b/src/maik/mxg.c @@ -0,0 +1,39 @@ +/* MK, 10.09.2010 + * + * This is a mex-wrapper for the computation of the function g of [Maischak]. */ + +#include +#include +#include "mex.h" +#include "G.h" + +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + const char* functionName = mexFunctionName(); + char errorMessage[255]; + char message[255]; + double y,x,a,k,p; + double* rval = NULL; + double* tmp = NULL; + int i; + + if (nlhs != 1 || nrhs != 5) { + sprintf(errorMessage, + "Error in %s: Please use exactly 1 output parameter and 5 " + "input parameters.", functionName); + mexErrMsgTxt(errorMessage); + return; + } + + y = mxGetScalar(prhs[0]); + x = mxGetScalar(prhs[1]); + a = mxGetScalar(prhs[2]); + k = mxGetScalar(prhs[3]); + p = mxGetScalar(prhs[4]); + + rval = compute_g(y,x,a,k,p); + plhs[0] = mxCreateDoubleMatrix(k+1,1,mxREAL); + tmp = mxGetPr(plhs[0]); + + for (i=0; i<=k; i++) + tmp[i] = rval[i]; +} diff --git a/src/maik/mxslpAD.c b/src/maik/mxslpAD.c new file mode 100644 index 0000000..5026437 --- /dev/null +++ b/src/maik/mxslpAD.c @@ -0,0 +1,47 @@ +/* MK, 10.09.2010 + * + * This is a mex-wrapper for the computation of the antiderivative of the + * simple-layer potential */ + +#include +#include +#include "mex.h" +#include "G.h" + +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + const char* functionName = mexFunctionName(); + char errorMessage[255]; + double y1,y2,x1,x2,a,k,ell; + double** rval =NULL; + double* tmp =NULL; + int i,j; + + if (nlhs != 1 || nrhs != 7) { + sprintf(errorMessage, + "Error in %s: Please use exactly 1 output parameter and 7 " + "input parameters.", functionName); + mexErrMsgTxt(errorMessage); + return; + } + + y1 = mxGetScalar(prhs[0]); + y2 = mxGetScalar(prhs[1]); + x1 = mxGetScalar(prhs[2]); + x2 = mxGetScalar(prhs[3]); + a = mxGetScalar(prhs[4]); + k = mxGetScalar(prhs[5]); + ell = mxGetScalar(prhs[6]); + + rval = slpAD(y1,y2,x1,x2,a,k,ell); + /* plhs[0] = mxCreateDoubleScalar((double) rval);*/ + plhs[0] = mxCreateDoubleMatrix(k+1,ell+1,mxREAL); + tmp = mxGetPr(plhs[0]); + + for ( i=0;i<=k; i++ ){ + for ( j=0; j<=ell; j++ ){ + tmp[j*((int)k+1)+i] = rval[i][j]; + } + } + + free(rval); +} diff --git a/src/maik/mxslpADLO.c b/src/maik/mxslpADLO.c new file mode 100644 index 0000000..76886c9 --- /dev/null +++ b/src/maik/mxslpADLO.c @@ -0,0 +1,33 @@ +/* MK, 10.09.2010 + * + * This is a mex-wrapper for the computation of the antiderivative of the + * simple-layer potential in the lowest-order case */ + +#include +#include +#include "mex.h" +#include "G.h" + +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + const char* functionName = mexFunctionName(); + char errorMessage[255]; + double y1,y2,x1,x2,a; + double rval =0; + + if (nlhs != 1 || nrhs != 5) { + sprintf(errorMessage, + "Error in %s: Please use exactly 1 output parameter and 5 " + "input parameters.", functionName); + mexErrMsgTxt(errorMessage); + return; + } + + y1 = mxGetScalar(prhs[0]); + y2 = mxGetScalar(prhs[1]); + x1 = mxGetScalar(prhs[2]); + x2 = mxGetScalar(prhs[3]); + a = mxGetScalar(prhs[4]); + + rval = slpADLO(y1,y2,x1,x2,a); + plhs[0] = mxCreateDoubleScalar((double) rval); +} diff --git a/src/maik/slpAD_test_A_C.c b/src/maik/slpAD_test_A_C.c new file mode 100644 index 0000000..188348c --- /dev/null +++ b/src/maik/slpAD_test_A_C.c @@ -0,0 +1,36 @@ +/* + * ===================================================================================== + * + * Filename: slpAD_text_A_C.c + * + * Description: + * + * Version: 1.0 + * Created: 09/10/2010 12:29:04 PM + * Revision: none + * Compiler: gcc + * + * Author: YOUR NAME (), + * Company: + * + * ===================================================================================== + */ + +#include "G.h" +#include + +void main(void) { + double a; + + a = slpAD(1,1,5,5,4,0,0); + a = slpAD(1,1,5,5,4,0,0); + a = slpAD(1,1,5,5,4,0,0); + a = slpAD(1,1,5,5,4,0,0); + a = slpAD(1,1,5,5,4,0,0); + a = slpAD(1,1,5,5,4,0,0); + a = slpAD(1,1,5,5,4,0,0); + a = slpAD(1,1,5,5,4,0,0); + a = slpAD(1,1,5,5,4,0,0); + a = slpAD(1,1,5,5,4,0,0); + printf("blabla %f",a); +} diff --git a/src/maik/test_g.m b/src/maik/test_g.m new file mode 100644 index 0000000..d46791a --- /dev/null +++ b/src/maik/test_g.m @@ -0,0 +1,65 @@ +function [fails,passes,relative_errors] = test_g() + +% test_g: Performs several tests to verify the function which computes g, +% needed for the simple-layer potential. +% +% Usage: [fails,passes,relative_errors] = test_g() +% +% Description: +% +%Output: fails ............. number of failed testcases. +% passes ............ number of passed testcases. +% relative_errors ... errors compaired to the reference solution +% +% Authors: Michael Karkulik +% +% (c) Hilbert-Team 2010 + +passes=0; +fails=0; +relative_errors = zeros(10,1); +kmax = 4; % maximal polynomial degree + +% p = -1/2 +p = -1/2; +x = -1 + [0:0.3:0.9]'; +nx = size(x,1); +a = [-5:1:5]'; +na = size(a,1); + +a = reshape(repmat(a',nx,1),nx*na,1); +x = repmat(x,na,1); + +for i = 1:size(a,1) + for k = 0:kmax + f=@(y)(integrand(y,x(i),a(i),k,p)); + reference(i,k+1) = quad(f,0,1,1e4*eps); + end + s = mxg(1,x(i),a(i),kmax,p); + t = mxg(0,x(i),a(i),kmax,p); + value(i,:) = s-t; +end +reference +value +relative_errors = abs(reference-value)./abs(reference) + +for j=1:size(relative_errors,1) + for i=1:size(relative_errors,2) + if (relative_errors(j,i)>=1e6*eps) + fails = fails+1; + else + passes=passes+1; + end + end +end + +end +%------------------------------------------------------------------------------ + +function val = integrand(y,x,a,k,p) +val = zeros(length(y),1); +val = (y.^k).*( (y-x).*(y-x) + a.*a).^p; +%for i=1:length(y) +% val(i) = y(i).^k.*( (y(i)-x).*(y(i)-x) + a.*a ).^p; +%end +end diff --git a/src/maik/test_slpAD.m b/src/maik/test_slpAD.m new file mode 100644 index 0000000..ed90b3b --- /dev/null +++ b/src/maik/test_slpAD.m @@ -0,0 +1,26 @@ +function [fails,passes,relative_errors] = test_slpAD() + +% test_slpAD: Performs tests for the computation of the function slpAD, which is +% used for the computation of the discrete simple-layer potential + +x1 = 2; +x2 = 2; +a = 0; +k=0; +ell=0; +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) +f = @(s,t) integrand(s,t,x1,x2,a,k,ell); +reference = dblquad(f,0,1,0,1,1e4*eps) + + + +end +%------------------------------------------------------------------------------ +function val = integrand(y1,y2,x1,x2,lambda,k,ell) +val = zeros(length(y1),length(y2)); +for i=1:length(y1) + for j=1:length(y2) + 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); + end +end +end diff --git a/src/maik/test_slpADLO.m b/src/maik/test_slpADLO.m new file mode 100644 index 0000000..de50dc1 --- /dev/null +++ b/src/maik/test_slpADLO.m @@ -0,0 +1,29 @@ +function [fails,passes,relative_errors] = test_slpADLO() + +% test_slpAD: Performs tests for the computation of the function slpADLO, which is +% used for the computation of the lowest +% order discrete simple-layer potential. + +x1 = 2; +x2 = 2; +a = 0; +A = mxslpADLO(0,0,x1,x2,a); +B = mxslpADLO(0,1,x1,x2,a); +C = mxslpADLO(1,0,x1,x2,a); +D = mxslpADLO(1,1,x1,x2,a); +value = A+D-B-C +f = @(s,t) integrand(s,t,x1,x2,a); +reference = dblquad(f,0,1,0,1,1e4*eps) + + + +end +%------------------------------------------------------------------------------ +function val = integrand(y1,y2,x1,x2,lambda) +val = zeros(length(y1),length(y2)); +for i=1:length(y1) + for j=1:length(y2) + val(i,j) = ( (x1-y1(i))*(x1-y1(i)) + (x2-y2(j))*(x2-y2(j)) + lambda*lambda )^(-1/2); + end +end +end -- 2.47.3