From: treecity Date: Wed, 11 May 2011 16:47:54 +0000 (+0000) Subject: git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@22 26120e32-c555-405d... X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=7b568236e2e40a08f9e9c5f90f3138c2a20a4f1d;p=bacc.git git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@22 26120e32-c555-405d-b3e1-1f783fb42516 --- diff --git a/src/MK/G.c b/src/MK/G.c deleted file mode 100644 index e0b6699..0000000 --- a/src/MK/G.c +++ /dev/null @@ -1,215 +0,0 @@ -#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/MK/G.h b/src/MK/G.h deleted file mode 100644 index aec9e09..0000000 --- a/src/MK/G.h +++ /dev/null @@ -1,23 +0,0 @@ -#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/MK/Makefile b/src/MK/Makefile deleted file mode 100644 index cee0636..0000000 --- a/src/MK/Makefile +++ /dev/null @@ -1,10 +0,0 @@ -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/MK/g.c b/src/MK/g.c deleted file mode 100644 index eb92f9c..0000000 --- a/src/MK/g.c +++ /dev/null @@ -1,70 +0,0 @@ -#include "g.h" - -/* MK: fast versions of the functios compute_g possible? */ - -double compute_g0(double y, double x, double a, double p) { - int sp = (int) 2*(p-EPS); // MK (p-EPS) instead of (p+EPS) - //printf("\n compute_g0, p = %lf, sp = %d\n",p,sp); - assert(p == 0 - || (p == -0.5) - || (p == -1) - || (p == -1.5) - || (fabs(a) <= EPS) ); - if (fabs(a) <= EPS){ - // printf("\n a < eps\n"); - switch(sp) { - case 0: - return y-x; - case -1: - return log(fabs(x-y))*(y-x)/fabs(y-x); - case -2: - return - (y-x)/fabs(y-x)/fabs(y-x); - case -3: - return -0.5*(y-x)/fabs(y-x)/fabs(y-x)/fabs(y-x); - } - } else { - // printf("\n a > eps\n"); - switch (sp) { - case 0: - return y-x; - case -1: - return asinh((y-x)/fabs(a)); - case -2: - return atan((y-x)/fabs(a)); - case -3: - return (y-x)*pow((x*x+y*y+a*a-2*x*y), -0.5)/(a*a); - default: - ERROR_MSG("p must be either 0, -1/2, -1 or -3/2."); - return NAN; - } - } -} - -// MK: the cases p=-1, p=-3/2 do not work!!! (division by zero after some polynomial degree) -double* compute_g(double y, double x, double a, int k, double p) { - int i = 0; - double* gk = (double*) malloc(sizeof(double) * (k+1)); - double factor = x*x+y*y+a*a-2*x*y; - double pot = pow(factor, p+1); - double ypow = y; - - - assert(p == 0 || p == -0.5 || p == -1 || p == -1.5); - - gk[0] = compute_g0(y,x,a,p); - - if (k >= 1) - gk[1] = (pot + 2.*x*(1.+p)*gk[0]) / (((double) 1)+2.*p+1.); - - for (i = 2; i <= k; ++i) { - gk[i] = (y*pot + 2.*x*(((double) i)+p)*gk[i-1] - (i-1)*(a*a+x*x)*gk[i-2]) - / (((double) i) + 2.*p + 1); - ypow *= y; - } - - return gk; -} - -/* History: - * MK 10.09.10 -> corrections, now the function seems to compute the right thing */ - diff --git a/src/MK/g.h b/src/MK/g.h deleted file mode 100644 index c950a56..0000000 --- a/src/MK/g.h +++ /dev/null @@ -1,22 +0,0 @@ -#ifndef _G_H_GUARD_ -#define _G_H_GUARD_ - -#include -#include -#include -#include - -#define EPS 1e-12 - -#if __STDC_VERSION__ >= 199901L -# define ERROR_MSG(msg) fprintf(stderr, "Error in %s at %s:%d:\n %s\n",\ - __func__, __FILE__, __LINE__, msg); -#else -# define ERROR_MSG(msg) fprintf(stderr, "Error at %s:%d:\n %s\n",\ - __FILE__, __LINE__, msg); -#endif - -double compute_g0(double y, double x, double a, double p); -double* compute_g(double y, double x, double a, int k, double p); - -#endif diff --git a/src/MK/mxg.c b/src/MK/mxg.c deleted file mode 100644 index c104919..0000000 --- a/src/MK/mxg.c +++ /dev/null @@ -1,39 +0,0 @@ -/* 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/MK/mxslpAD.c b/src/MK/mxslpAD.c deleted file mode 100644 index 5026437..0000000 --- a/src/MK/mxslpAD.c +++ /dev/null @@ -1,47 +0,0 @@ -/* 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/MK/mxslpADLO.c b/src/MK/mxslpADLO.c deleted file mode 100644 index 76886c9..0000000 --- a/src/MK/mxslpADLO.c +++ /dev/null @@ -1,33 +0,0 @@ -/* 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/MK/slpAD_test_A_C.c b/src/MK/slpAD_test_A_C.c deleted file mode 100644 index 188348c..0000000 --- a/src/MK/slpAD_test_A_C.c +++ /dev/null @@ -1,36 +0,0 @@ -/* - * ===================================================================================== - * - * 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/MK/test_g.m b/src/MK/test_g.m deleted file mode 100644 index d46791a..0000000 --- a/src/MK/test_g.m +++ /dev/null @@ -1,65 +0,0 @@ -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/MK/test_slpAD.m b/src/MK/test_slpAD.m deleted file mode 100644 index ed90b3b..0000000 --- a/src/MK/test_slpAD.m +++ /dev/null @@ -1,26 +0,0 @@ -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/MK/test_slpADLO.m b/src/MK/test_slpADLO.m deleted file mode 100644 index de50dc1..0000000 --- a/src/MK/test_slpADLO.m +++ /dev/null @@ -1,29 +0,0 @@ -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