From: treecity Date: Wed, 11 May 2011 05:50:13 +0000 (+0000) Subject: [bem3d] neuste Version eingescheckt X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=2641b4b7beebcf59a20eae9ff1d35ee7e620139e;p=bacc.git [bem3d] neuste Version eingescheckt [src] test_g0 -> g0 einzeln getestet (gegen Karkulik und Quadratur) git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@20 26120e32-c555-405d-b3e1-1f783fb42516 --- diff --git a/src/MK/G.c b/src/MK/G.c new file mode 100644 index 0000000..e0b6699 --- /dev/null +++ b/src/MK/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/MK/G.h b/src/MK/G.h new file mode 100644 index 0000000..aec9e09 --- /dev/null +++ b/src/MK/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/MK/Makefile b/src/MK/Makefile new file mode 100644 index 0000000..cee0636 --- /dev/null +++ b/src/MK/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/MK/g.c b/src/MK/g.c new file mode 100644 index 0000000..eb92f9c --- /dev/null +++ b/src/MK/g.c @@ -0,0 +1,70 @@ +#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 new file mode 100644 index 0000000..c950a56 --- /dev/null +++ b/src/MK/g.h @@ -0,0 +1,22 @@ +#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 new file mode 100644 index 0000000..c104919 --- /dev/null +++ b/src/MK/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/MK/mxslpAD.c b/src/MK/mxslpAD.c new file mode 100644 index 0000000..5026437 --- /dev/null +++ b/src/MK/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/MK/mxslpADLO.c b/src/MK/mxslpADLO.c new file mode 100644 index 0000000..76886c9 --- /dev/null +++ b/src/MK/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/MK/slpAD_test_A_C.c b/src/MK/slpAD_test_A_C.c new file mode 100644 index 0000000..188348c --- /dev/null +++ b/src/MK/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/MK/test_g.m b/src/MK/test_g.m new file mode 100644 index 0000000..d46791a --- /dev/null +++ b/src/MK/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/MK/test_slpAD.m b/src/MK/test_slpAD.m new file mode 100644 index 0000000..ed90b3b --- /dev/null +++ b/src/MK/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/MK/test_slpADLO.m b/src/MK/test_slpADLO.m new file mode 100644 index 0000000..de50dc1 --- /dev/null +++ b/src/MK/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 diff --git a/src/SLPrecangle.cpp b/src/SLPrecangle.cpp new file mode 100644 index 0000000..5a414b9 --- /dev/null +++ b/src/SLPrecangle.cpp @@ -0,0 +1,358 @@ +#include +#include +#include +#include +#include "mex.h" + +#include "SLPrecangle.hpp" + +#define my_PI 3.141592653589793 +#define EPS 0.00001 + +using namespace std; + +int sign(double); +double arsinh(double); + +// sol = quadInt((F_par/F_ort),s1,s2,k1,k2,t1,t2,l1,l2,d1,d2,d3); +double quadInt( + double(*)(double, double, double, double, double, double, double), + double, double, double, double, double, double, double, double, double, + double, double); +/* +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { + //sicherheitsabfragen zu Datengroessen + if (nrhs != 2) + mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx3))"); + if (nlhs > 1) + mexErrMsgTxt("has only one output argument"); + + int cm = mxGetM(prhs[0]); + int cn = mxGetN(prhs[0]); + + if (cn != 3) + mexErrMsgTxt("expected coordinates (Nx3)"); + int em = mxGetM(prhs[1]); + int en = mxGetN(prhs[1]); + if (en != 3) + mexErrMsgTxt("expected elements (Mx3)"); + //Vorbereitung der Daten + + plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL); + double * A = mxGetPr(plhs[0]); + double * C = mxGetPr(prhs[0]); + double * E = mxGetPr(prhs[1]); + + double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * xn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + + double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + double * yn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + + double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); + + double tmp; + + int rx, ry; + + //Ausrechnen + for (int j = 0; j < em; ++j) { + x1[0] = C[(int) E[j] - 1]; + x1[1] = C[cm + (int) E[j] - 1]; + x1[2] = C[2 * cm + (int) E[j] - 1]; + + x2[0] = C[(int) E[em + j] - 1]; + x2[1] = C[cm + (int) E[em + j] - 1]; + x2[2] = C[2 * cm + (int) E[em + j] - 1]; + + x3[0] = C[(int) E[2 * em + j] - 1]; + x3[1] = C[cm + (int) E[2 * em + j] - 1]; + x3[2] = C[2 * cm + (int) E[2 * em + j] - 1]; + + xn[0] = (x2[1] - x1[1]) * (x3[2] - x1[2]) - (x2[2] - x1[2]) * (x3[1] + - x1[1]); + xn[1] = (x2[2] - x1[2]) * (x3[0] - x1[0]) - (x2[0] - x1[0]) * (x3[2] + - x1[2]); + xn[2] = (x2[0] - x1[0]) * (x3[1] - x1[1]) - (x2[1] - x1[1]) * (x3[0] + - x1[0]); + + //if(xn[0]*xn[0]+xn[1]*xn[1] 0 ? 1 : (x < 0 ? -1 : 0); +} + +double inline arsinh(double x) { + return log(x + sqrt(x * x + 1)); +} + +//y-x muss != 0 sein +double g0(double p, double y, double x, double l) { + //printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l); + + double sol = 0; + + if (l != 0) { + if (p == 0.5) { + sol = (y - x) / 2 * sqrt((y - x) * (y - x) + l * l) + l * l / 2 + * asinh((y - x) / fabs(l)); + // printf("%.2f |",sol); + } else if (p == 0) + sol = y - x; + else if (p == -0.5) + sol = asinh((y - x) / fabs(l)); + else if (p == -1) + sol = atan((y - x) / fabs(l)); + else if (p == -1.5) + sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l)); + else + sol = (y - x) * pow((y - x) * (y - x) + l * l, p) + 2 * p * l * l + * g0(p - 1, y, x, l) / (2 * p + 1); + } else { + if (p == -0.5) + sol = sign(y - x) * log(fabs(y - x)); + else + sol = (y - x) * pow(fabs(y - x), 2 * p) / (2 * p + 1); + } + + return sol; +} + +double G00(double p, double y1, double y2, double x1, double x2, double l) { + // printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l); + + double sol = 0; + if (p == -1.5) { + if (l == 0) { + sol = -sqrt((y1 - x1) * (y1 - x1) + (y2 - x2) * (y2 - x2)) / ((y1 + - x1) * (y2 - x2)); + } else { + sol = sign((y1 - x1) * (y2 - x2)) / (2 * fabs(l)) * acos( + -2 * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2) / (((y1 + - x1) * (y1 - x1) + l * l) * ((y2 - x2) * (y2 - x2) + + l * l)) + 1); + } + + } else if (p == -0.5) { + if (l != 0) + sol = 2 * p * l * l * G00(p - 1, y1, y2, x1, x2, l); + if ((y1 - x1) != 0) + sol += (y1 - x1) * g0(p, y2, x2, + sqrt((y1 - x1) * (y1 - x1) + l * l)); + if ((y2 - x2) != 0) + sol += (y2 - x2) * g0(p, y1, x1, + sqrt((y2 - x2) * (y2 - x2) + l * l)); + sol /= 2 * p + 2; + } else { + mexErrMsgTxt("no case for p defined"); + } + + return sol; +} + +double F_par(double x1, double x2, double y1, double y2, double d1, double d2, + double d3) { + + // printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3); + double sol = (x1 - y1 - d1) * (x2 - y2 - d2); + + if (sol != 0) + sol *= slpADLO(x1, x2, y1 + d1, y2 + d2, d3); + + if ((x1 - y1 - d1) != 0) + sol -= (x1 - y1 - d1) * compute_g0(0.5, x1, y1 + d1, + sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3)); + + if ((x2 - y2 - d2) != 0) + sol -= (x2 - y2 - d2) * compute_g0(0.5, x2, y2 + d2, + sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3)); + + double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2 + - d2) + d3 * d3); + sol += 1. / 3 * hlp * sqrt(hlp); + return sol; +} + +double inline quadInt( + double(*f)(double, double, double, double, double, double, double), + double s1, double s2, double k1, double k2, double t1, double t2, + double l1, double l2, double d1, double d2, double d3) { + + return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f( + k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1, + s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2, + t1, l2, d1, d2, d3) - f(k1, s2, t1, t2, d1, d2, d3) - f(s1, k2, l1, + l2, d1, d2, d3) + f(s1, k2, l1, t2, d1, d2, d3) + f(s1, k2, t1, l2, + d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1, + d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2, + d3) + f(s1, s2, t1, t2, d1, d2, d3); +} + +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(-0.5, y2, x2, tmp); + } + tmp = y2 - x2; + if (fabs(tmp) >= EPS * y2) { + tmp = sqrt(y2 * y2 + x2 * x2 + a * a - 2 * y2 * x2); + gK = compute_g0(-0.5, y1, x1, tmp); + } + 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; +} + +double compute_g0(double p, double y, double x, double a) { + 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: + printf("p must be either 0, -1/2, -1 or -3/2."); + return NAN; + } + } +} + +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(-0.5, + x1, yd1, tmp3) - tmp2 * compute_g0(-0.5, x2, yd2, tmp4) + tmp5 / 3.; + return rval; +} diff --git a/src/SLPrecangle.hpp b/src/SLPrecangle.hpp new file mode 100644 index 0000000..6f80150 --- /dev/null +++ b/src/SLPrecangle.hpp @@ -0,0 +1,21 @@ +#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 diff --git a/src/bem3d/ROADMAP b/src/bem3d/ROADMAP index f32dec8..8c8b144 100644 --- a/src/bem3d/ROADMAP +++ b/src/bem3d/ROADMAP @@ -22,7 +22,7 @@ | 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 | | | @@ -47,16 +47,11 @@ | 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. ******************************************************************************* diff --git a/src/bem3d/src/boundary_mesh/HalfEdge.hpp b/src/bem3d/src/boundary_mesh/HalfEdge.hpp index 6d3193b..9f1f222 100644 --- a/src/bem3d/src/boundary_mesh/HalfEdge.hpp +++ b/src/bem3d/src/boundary_mesh/HalfEdge.hpp @@ -44,6 +44,8 @@ namespace boundary_mesh { 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; diff --git a/src/bem3d/src/boundary_mesh/HalfEdge_tmpl.cpp b/src/bem3d/src/boundary_mesh/HalfEdge_tmpl.cpp index f62a197..8fb7b29 100644 --- a/src/bem3d/src/boundary_mesh/HalfEdge_tmpl.cpp +++ b/src/bem3d/src/boundary_mesh/HalfEdge_tmpl.cpp @@ -1,6 +1,32 @@ #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 { diff --git a/src/bem3d/t/generic/GenericVector.cpp b/src/bem3d/t/generic/GenericVector.cpp new file mode 100644 index 0000000..77f0ea4 --- /dev/null +++ b/src/bem3d/t/generic/GenericVector.cpp @@ -0,0 +1,10 @@ +#include +#include + +#define BOOST_TEST_MODULE generic::GenericVector test +#define BOOST_TEST_DYN_LINK + +#include + +#include "GenericVector/01-basics.cpp" + diff --git a/src/bem3d/t/generic/GenericVector/01-basics.cpp b/src/bem3d/t/generic/GenericVector/01-basics.cpp new file mode 100644 index 0000000..bc2bae4 --- /dev/null +++ b/src/bem3d/t/generic/GenericVector/01-basics.cpp @@ -0,0 +1,26 @@ +#include + +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 ); +} + diff --git a/src/mex_g0.cpp b/src/mex_g0.cpp new file mode 100644 index 0000000..b31afe4 --- /dev/null +++ b/src/mex_g0.cpp @@ -0,0 +1,36 @@ +#include +#include +#include +#include "mex.h" +#include + +#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])); + + +} + + + diff --git a/src/test_g0.m b/src/test_g0.m new file mode 100644 index 0000000..b08aaf0 --- /dev/null +++ b/src/test_g0.m @@ -0,0 +1,16 @@ +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 +