diff --git a/matlab/fmm2d.c b/matlab/fmm2d.c index 22eb84d..6b2b317 100644 --- a/matlab/fmm2d.c +++ b/matlab/fmm2d.c @@ -2,7 +2,7 @@ /* Automatically generated by mwrap */ /* --------------------------------------------------- */ -/* Code generated by mwrap 1.0 */ +/* Code generated by mwrap 1.1 */ /* Copyright statement for mwrap: @@ -1104,6 +1104,9 @@ mxWrapReturnZDef_single (mxWrapReturn_single_dcomplex, dcomplex, #define MWF77_c2d_directch C2D_DIRECTCH #define MWF77_c2d_directdh C2D_DIRECTDH #define MWF77_c2d_directcdh C2D_DIRECTCDH +#define MWF77_stfmm2d STFMM2D +#define MWF77_st2ddirectstokg ST2DDIRECTSTOKG +#define MWF77_st2ddirectstokstrsg ST2DDIRECTSTOKSTRSG #elif defined(MWF77_UNDERSCORE1) #define MWF77_hndiv2d hndiv2d_ #define MWF77_hfmm2d_ndiv hfmm2d_ndiv_ @@ -1146,6 +1149,9 @@ mxWrapReturnZDef_single (mxWrapReturn_single_dcomplex, dcomplex, #define MWF77_c2d_directch c2d_directch_ #define MWF77_c2d_directdh c2d_directdh_ #define MWF77_c2d_directcdh c2d_directcdh_ +#define MWF77_stfmm2d stfmm2d_ +#define MWF77_st2ddirectstokg st2ddirectstokg_ +#define MWF77_st2ddirectstokstrsg st2ddirectstokstrsg_ #elif defined(MWF77_UNDERSCORE0) #define MWF77_hndiv2d hndiv2d #define MWF77_hfmm2d_ndiv hfmm2d_ndiv @@ -1188,6 +1194,9 @@ mxWrapReturnZDef_single (mxWrapReturn_single_dcomplex, dcomplex, #define MWF77_c2d_directch c2d_directch #define MWF77_c2d_directdh c2d_directdh #define MWF77_c2d_directcdh c2d_directcdh +#define MWF77_stfmm2d stfmm2d +#define MWF77_st2ddirectstokg st2ddirectstokg +#define MWF77_st2ddirectstokstrsg st2ddirectstokstrsg #else /* f2c convention */ #define MWF77_hndiv2d hndiv2d_ #define MWF77_hfmm2d_ndiv hfmm2d_ndiv__ @@ -1230,6 +1239,9 @@ mxWrapReturnZDef_single (mxWrapReturn_single_dcomplex, dcomplex, #define MWF77_c2d_directch c2d_directch__ #define MWF77_c2d_directdh c2d_directdh__ #define MWF77_c2d_directcdh c2d_directcdh__ +#define MWF77_stfmm2d stfmm2d_ +#define MWF77_st2ddirectstokg st2ddirectstokg_ +#define MWF77_st2ddirectstokstrsg st2ddirectstokstrsg_ #endif #ifdef __cplusplus @@ -1281,6 +1293,9 @@ MWF77_RETURN MWF77_c2d_directcdg(int*, double*, int*, dcomplex*, dcomplex*, doub MWF77_RETURN MWF77_c2d_directch(int*, double*, int*, dcomplex*, double*, int*, dcomplex*, dcomplex*, dcomplex*, double*); MWF77_RETURN MWF77_c2d_directdh(int*, double*, int*, dcomplex*, double*, int*, dcomplex*, dcomplex*, dcomplex*, double*); MWF77_RETURN MWF77_c2d_directcdh(int*, double*, int*, dcomplex*, dcomplex*, double*, int*, dcomplex*, dcomplex*, dcomplex*, double*); +MWF77_RETURN MWF77_stfmm2d(int*, double*, int*, double*, int*, double*, int*, double*, double*, int*, double*, double*, double*, int*, double*, int*, double*, double*, double*, int*); +MWF77_RETURN MWF77_st2ddirectstokg(int*, double*, double*, int*, double*, int*, double*, double*, double*, double*); +MWF77_RETURN MWF77_st2ddirectstokstrsg(int*, double*, int*, double*, int*, double*, double*, int*, double*, int*, double*, double*, double*, double*); #ifdef __cplusplus } /* end extern C */ @@ -11089,6 +11104,872 @@ void mexStub44(int nlhs, mxArray* plhs[], mexErrMsgTxt(mw_err_txt_); } +/* ---- fmm2d.mw: 1647 ---- + * stfmm2d(int[1] nd, double[1] eps, int[1] ns, double[2, ns] sources, int[1] ifstoklet, double[nd2, ns_stok] stoklet, int[1] ifstrslet, double[nd2, ns_strs] strslet, double[nd2, ns_strs] strsvec, int[1] ifppreg, inout double[nd2, ns_pot] pot, inout double[nd, ns_pre] pre, inout double[nd4, ns_grad] grad, int[1] nt, double[2, nt] targ, int[1] ifppregtarg, inout double[nd2, nt_pot] pottarg, inout double[nd, nt_pre] pretarg, inout double[nd4, nt_grad] gradtarg, inout int[1] ier); + */ +static const char* stubids45_ = "stfmm2d(i int[x], i double[x], i int[x], i double[xx], i int[x], i double[xx], i int[x], i double[xx], i double[xx], i int[x], io double[xx], io double[xx], io double[xx], i int[x], i double[xx], i int[x], io double[xx], io double[xx], io double[xx], io int[x])"; + +void mexStub45(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + const char* mw_err_txt_ = 0; + int* in0_ =0; /* nd */ + double* in1_ =0; /* eps */ + int* in2_ =0; /* ns */ + double* in3_ =0; /* sources */ + int* in4_ =0; /* ifstoklet */ + double* in5_ =0; /* stoklet */ + int* in6_ =0; /* ifstrslet */ + double* in7_ =0; /* strslet */ + double* in8_ =0; /* strsvec */ + int* in9_ =0; /* ifppreg */ + double* in10_ =0; /* pot */ + double* in11_ =0; /* pre */ + double* in12_ =0; /* grad */ + int* in13_ =0; /* nt */ + double* in14_ =0; /* targ */ + int* in15_ =0; /* ifppregtarg */ + double* in16_ =0; /* pottarg */ + double* in17_ =0; /* pretarg */ + double* in18_ =0; /* gradtarg */ + int* in19_ =0; /* ier */ + mwSize dim20_; /* 1 */ + mwSize dim21_; /* 1 */ + mwSize dim22_; /* 1 */ + mwSize dim23_; /* 2 */ + mwSize dim24_; /* ns */ + mwSize dim25_; /* 1 */ + mwSize dim26_; /* nd2 */ + mwSize dim27_; /* ns_stok */ + mwSize dim28_; /* 1 */ + mwSize dim29_; /* nd2 */ + mwSize dim30_; /* ns_strs */ + mwSize dim31_; /* nd2 */ + mwSize dim32_; /* ns_strs */ + mwSize dim33_; /* 1 */ + mwSize dim34_; /* nd2 */ + mwSize dim35_; /* ns_pot */ + mwSize dim36_; /* nd */ + mwSize dim37_; /* ns_pre */ + mwSize dim38_; /* nd4 */ + mwSize dim39_; /* ns_grad */ + mwSize dim40_; /* 1 */ + mwSize dim41_; /* 2 */ + mwSize dim42_; /* nt */ + mwSize dim43_; /* 1 */ + mwSize dim44_; /* nd2 */ + mwSize dim45_; /* nt_pot */ + mwSize dim46_; /* nd */ + mwSize dim47_; /* nt_pre */ + mwSize dim48_; /* nd4 */ + mwSize dim49_; /* nt_grad */ + mwSize dim50_; /* 1 */ + + dim20_ = (mwSize) mxWrapGetScalar(prhs[20], &mw_err_txt_); + dim21_ = (mwSize) mxWrapGetScalar(prhs[21], &mw_err_txt_); + dim22_ = (mwSize) mxWrapGetScalar(prhs[22], &mw_err_txt_); + dim23_ = (mwSize) mxWrapGetScalar(prhs[23], &mw_err_txt_); + dim24_ = (mwSize) mxWrapGetScalar(prhs[24], &mw_err_txt_); + dim25_ = (mwSize) mxWrapGetScalar(prhs[25], &mw_err_txt_); + dim26_ = (mwSize) mxWrapGetScalar(prhs[26], &mw_err_txt_); + dim27_ = (mwSize) mxWrapGetScalar(prhs[27], &mw_err_txt_); + dim28_ = (mwSize) mxWrapGetScalar(prhs[28], &mw_err_txt_); + dim29_ = (mwSize) mxWrapGetScalar(prhs[29], &mw_err_txt_); + dim30_ = (mwSize) mxWrapGetScalar(prhs[30], &mw_err_txt_); + dim31_ = (mwSize) mxWrapGetScalar(prhs[31], &mw_err_txt_); + dim32_ = (mwSize) mxWrapGetScalar(prhs[32], &mw_err_txt_); + dim33_ = (mwSize) mxWrapGetScalar(prhs[33], &mw_err_txt_); + dim34_ = (mwSize) mxWrapGetScalar(prhs[34], &mw_err_txt_); + dim35_ = (mwSize) mxWrapGetScalar(prhs[35], &mw_err_txt_); + dim36_ = (mwSize) mxWrapGetScalar(prhs[36], &mw_err_txt_); + dim37_ = (mwSize) mxWrapGetScalar(prhs[37], &mw_err_txt_); + dim38_ = (mwSize) mxWrapGetScalar(prhs[38], &mw_err_txt_); + dim39_ = (mwSize) mxWrapGetScalar(prhs[39], &mw_err_txt_); + dim40_ = (mwSize) mxWrapGetScalar(prhs[40], &mw_err_txt_); + dim41_ = (mwSize) mxWrapGetScalar(prhs[41], &mw_err_txt_); + dim42_ = (mwSize) mxWrapGetScalar(prhs[42], &mw_err_txt_); + dim43_ = (mwSize) mxWrapGetScalar(prhs[43], &mw_err_txt_); + dim44_ = (mwSize) mxWrapGetScalar(prhs[44], &mw_err_txt_); + dim45_ = (mwSize) mxWrapGetScalar(prhs[45], &mw_err_txt_); + dim46_ = (mwSize) mxWrapGetScalar(prhs[46], &mw_err_txt_); + dim47_ = (mwSize) mxWrapGetScalar(prhs[47], &mw_err_txt_); + dim48_ = (mwSize) mxWrapGetScalar(prhs[48], &mw_err_txt_); + dim49_ = (mwSize) mxWrapGetScalar(prhs[49], &mw_err_txt_); + dim50_ = (mwSize) mxWrapGetScalar(prhs[50], &mw_err_txt_); + + if (mxGetM(prhs[0])*mxGetN(prhs[0]) != dim20_) { + mw_err_txt_ = "Bad argument size: nd"; goto mw_err_label; + } + + if (mxGetM(prhs[1])*mxGetN(prhs[1]) != dim21_) { + mw_err_txt_ = "Bad argument size: eps"; goto mw_err_label; + } + + if (mxGetM(prhs[2])*mxGetN(prhs[2]) != dim22_) { + mw_err_txt_ = "Bad argument size: ns"; goto mw_err_label; + } + + if (mxGetM(prhs[3]) != dim23_ || + mxGetN(prhs[3]) != dim24_) { + mw_err_txt_ = "Bad argument size: sources"; + goto mw_err_label; + } + + if (mxGetM(prhs[4])*mxGetN(prhs[4]) != dim25_) { + mw_err_txt_ = "Bad argument size: ifstoklet"; goto mw_err_label; + } + + if (mxGetM(prhs[5]) != dim26_ || + mxGetN(prhs[5]) != dim27_) { + mw_err_txt_ = "Bad argument size: stoklet"; + goto mw_err_label; + } + + if (mxGetM(prhs[6])*mxGetN(prhs[6]) != dim28_) { + mw_err_txt_ = "Bad argument size: ifstrslet"; goto mw_err_label; + } + + if (mxGetM(prhs[7]) != dim29_ || + mxGetN(prhs[7]) != dim30_) { + mw_err_txt_ = "Bad argument size: strslet"; + goto mw_err_label; + } + + if (mxGetM(prhs[8]) != dim31_ || + mxGetN(prhs[8]) != dim32_) { + mw_err_txt_ = "Bad argument size: strsvec"; + goto mw_err_label; + } + + if (mxGetM(prhs[9])*mxGetN(prhs[9]) != dim33_) { + mw_err_txt_ = "Bad argument size: ifppreg"; goto mw_err_label; + } + + if (mxGetM(prhs[10]) != dim34_ || + mxGetN(prhs[10]) != dim35_) { + mw_err_txt_ = "Bad argument size: pot"; + goto mw_err_label; + } + + if (mxGetM(prhs[11]) != dim36_ || + mxGetN(prhs[11]) != dim37_) { + mw_err_txt_ = "Bad argument size: pre"; + goto mw_err_label; + } + + if (mxGetM(prhs[12]) != dim38_ || + mxGetN(prhs[12]) != dim39_) { + mw_err_txt_ = "Bad argument size: grad"; + goto mw_err_label; + } + + if (mxGetM(prhs[13])*mxGetN(prhs[13]) != dim40_) { + mw_err_txt_ = "Bad argument size: nt"; goto mw_err_label; + } + + if (mxGetM(prhs[14]) != dim41_ || + mxGetN(prhs[14]) != dim42_) { + mw_err_txt_ = "Bad argument size: targ"; + goto mw_err_label; + } + + if (mxGetM(prhs[15])*mxGetN(prhs[15]) != dim43_) { + mw_err_txt_ = "Bad argument size: ifppregtarg"; goto mw_err_label; + } + + if (mxGetM(prhs[16]) != dim44_ || + mxGetN(prhs[16]) != dim45_) { + mw_err_txt_ = "Bad argument size: pottarg"; + goto mw_err_label; + } + + if (mxGetM(prhs[17]) != dim46_ || + mxGetN(prhs[17]) != dim47_) { + mw_err_txt_ = "Bad argument size: pretarg"; + goto mw_err_label; + } + + if (mxGetM(prhs[18]) != dim48_ || + mxGetN(prhs[18]) != dim49_) { + mw_err_txt_ = "Bad argument size: gradtarg"; + goto mw_err_label; + } + + if (mxGetM(prhs[19])*mxGetN(prhs[19]) != dim50_) { + mw_err_txt_ = "Bad argument size: ier"; goto mw_err_label; + } + + if (mxGetM(prhs[0])*mxGetN(prhs[0]) != 0) { + in0_ = mxWrapGetArray_int(prhs[0], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in0_ = NULL; + if (mxGetM(prhs[1])*mxGetN(prhs[1]) != 0) { + if( mxGetClassID(prhs[1]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in1_ = mxGetDoubles(prhs[1]); +#else + in1_ = mxGetPr(prhs[1]); +#endif + } else + in1_ = NULL; + if (mxGetM(prhs[2])*mxGetN(prhs[2]) != 0) { + in2_ = mxWrapGetArray_int(prhs[2], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in2_ = NULL; + if (mxGetM(prhs[3])*mxGetN(prhs[3]) != 0) { + if( mxGetClassID(prhs[3]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in3_ = mxGetDoubles(prhs[3]); +#else + in3_ = mxGetPr(prhs[3]); +#endif + } else + in3_ = NULL; + if (mxGetM(prhs[4])*mxGetN(prhs[4]) != 0) { + in4_ = mxWrapGetArray_int(prhs[4], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in4_ = NULL; + if (mxGetM(prhs[5])*mxGetN(prhs[5]) != 0) { + if( mxGetClassID(prhs[5]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in5_ = mxGetDoubles(prhs[5]); +#else + in5_ = mxGetPr(prhs[5]); +#endif + } else + in5_ = NULL; + if (mxGetM(prhs[6])*mxGetN(prhs[6]) != 0) { + in6_ = mxWrapGetArray_int(prhs[6], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in6_ = NULL; + if (mxGetM(prhs[7])*mxGetN(prhs[7]) != 0) { + if( mxGetClassID(prhs[7]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in7_ = mxGetDoubles(prhs[7]); +#else + in7_ = mxGetPr(prhs[7]); +#endif + } else + in7_ = NULL; + if (mxGetM(prhs[8])*mxGetN(prhs[8]) != 0) { + if( mxGetClassID(prhs[8]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in8_ = mxGetDoubles(prhs[8]); +#else + in8_ = mxGetPr(prhs[8]); +#endif + } else + in8_ = NULL; + if (mxGetM(prhs[9])*mxGetN(prhs[9]) != 0) { + in9_ = mxWrapGetArray_int(prhs[9], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in9_ = NULL; + if (mxGetM(prhs[10])*mxGetN(prhs[10]) != 0) { + in10_ = mxWrapGetArray_double(prhs[10], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in10_ = NULL; + if (mxGetM(prhs[11])*mxGetN(prhs[11]) != 0) { + in11_ = mxWrapGetArray_double(prhs[11], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in11_ = NULL; + if (mxGetM(prhs[12])*mxGetN(prhs[12]) != 0) { + in12_ = mxWrapGetArray_double(prhs[12], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in12_ = NULL; + if (mxGetM(prhs[13])*mxGetN(prhs[13]) != 0) { + in13_ = mxWrapGetArray_int(prhs[13], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in13_ = NULL; + if (mxGetM(prhs[14])*mxGetN(prhs[14]) != 0) { + if( mxGetClassID(prhs[14]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in14_ = mxGetDoubles(prhs[14]); +#else + in14_ = mxGetPr(prhs[14]); +#endif + } else + in14_ = NULL; + if (mxGetM(prhs[15])*mxGetN(prhs[15]) != 0) { + in15_ = mxWrapGetArray_int(prhs[15], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in15_ = NULL; + if (mxGetM(prhs[16])*mxGetN(prhs[16]) != 0) { + in16_ = mxWrapGetArray_double(prhs[16], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in16_ = NULL; + if (mxGetM(prhs[17])*mxGetN(prhs[17]) != 0) { + in17_ = mxWrapGetArray_double(prhs[17], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in17_ = NULL; + if (mxGetM(prhs[18])*mxGetN(prhs[18]) != 0) { + in18_ = mxWrapGetArray_double(prhs[18], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in18_ = NULL; + if (mxGetM(prhs[19])*mxGetN(prhs[19]) != 0) { + in19_ = mxWrapGetArray_int(prhs[19], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in19_ = NULL; + if (mexprofrecord_) + mexprofrecord_[45]++; + MWF77_stfmm2d(in0_, in1_, in2_, in3_, in4_, in5_, in6_, in7_, in8_, in9_, in10_, in11_, in12_, in13_, in14_, in15_, in16_, in17_, in18_, in19_); + plhs[0] = mxCreateDoubleMatrix(dim34_, dim35_, mxREAL); + mxWrapCopy_double(plhs[0], in10_, dim34_*dim35_); + plhs[1] = mxCreateDoubleMatrix(dim36_, dim37_, mxREAL); + mxWrapCopy_double(plhs[1], in11_, dim36_*dim37_); + plhs[2] = mxCreateDoubleMatrix(dim38_, dim39_, mxREAL); + mxWrapCopy_double(plhs[2], in12_, dim38_*dim39_); + plhs[3] = mxCreateDoubleMatrix(dim44_, dim45_, mxREAL); + mxWrapCopy_double(plhs[3], in16_, dim44_*dim45_); + plhs[4] = mxCreateDoubleMatrix(dim46_, dim47_, mxREAL); + mxWrapCopy_double(plhs[4], in17_, dim46_*dim47_); + plhs[5] = mxCreateDoubleMatrix(dim48_, dim49_, mxREAL); + mxWrapCopy_double(plhs[5], in18_, dim48_*dim49_); + plhs[6] = mxCreateDoubleMatrix(dim50_, 1, mxREAL); + mxWrapCopy_int(plhs[6], in19_, dim50_); + +mw_err_label: + if (in0_) mxFree(in0_); + if (in2_) mxFree(in2_); + if (in4_) mxFree(in4_); + if (in6_) mxFree(in6_); + if (in9_) mxFree(in9_); + if (in10_) mxFree(in10_); + if (in11_) mxFree(in11_); + if (in12_) mxFree(in12_); + if (in13_) mxFree(in13_); + if (in15_) mxFree(in15_); + if (in16_) mxFree(in16_); + if (in17_) mxFree(in17_); + if (in18_) mxFree(in18_); + if (in19_) mxFree(in19_); + if (mw_err_txt_) + mexErrMsgTxt(mw_err_txt_); +} + +/* ---- fmm2d.mw: 1744 ---- + * st2ddirectstokg(int[1] nd, double[2, ns] sources, double[nd2, ns_stok] stoklet, int[1] ns, double[2, nt] targ, int[1] nt, inout double[nd2, nt] pottarg, inout double[nd, nt] pretarg, inout double[nd4, nt] gradtarg, double[1] thresh); + */ +static const char* stubids46_ = "st2ddirectstokg(i int[x], i double[xx], i double[xx], i int[x], i double[xx], i int[x], io double[xx], io double[xx], io double[xx], i double[x])"; + +void mexStub46(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + const char* mw_err_txt_ = 0; + int* in0_ =0; /* nd */ + double* in1_ =0; /* sources */ + double* in2_ =0; /* stoklet */ + int* in3_ =0; /* ns */ + double* in4_ =0; /* targ */ + int* in5_ =0; /* nt */ + double* in6_ =0; /* pottarg */ + double* in7_ =0; /* pretarg */ + double* in8_ =0; /* gradtarg */ + double* in9_ =0; /* thresh */ + mwSize dim10_; /* 1 */ + mwSize dim11_; /* 2 */ + mwSize dim12_; /* ns */ + mwSize dim13_; /* nd2 */ + mwSize dim14_; /* ns_stok */ + mwSize dim15_; /* 1 */ + mwSize dim16_; /* 2 */ + mwSize dim17_; /* nt */ + mwSize dim18_; /* 1 */ + mwSize dim19_; /* nd2 */ + mwSize dim20_; /* nt */ + mwSize dim21_; /* nd */ + mwSize dim22_; /* nt */ + mwSize dim23_; /* nd4 */ + mwSize dim24_; /* nt */ + mwSize dim25_; /* 1 */ + + dim10_ = (mwSize) mxWrapGetScalar(prhs[10], &mw_err_txt_); + dim11_ = (mwSize) mxWrapGetScalar(prhs[11], &mw_err_txt_); + dim12_ = (mwSize) mxWrapGetScalar(prhs[12], &mw_err_txt_); + dim13_ = (mwSize) mxWrapGetScalar(prhs[13], &mw_err_txt_); + dim14_ = (mwSize) mxWrapGetScalar(prhs[14], &mw_err_txt_); + dim15_ = (mwSize) mxWrapGetScalar(prhs[15], &mw_err_txt_); + dim16_ = (mwSize) mxWrapGetScalar(prhs[16], &mw_err_txt_); + dim17_ = (mwSize) mxWrapGetScalar(prhs[17], &mw_err_txt_); + dim18_ = (mwSize) mxWrapGetScalar(prhs[18], &mw_err_txt_); + dim19_ = (mwSize) mxWrapGetScalar(prhs[19], &mw_err_txt_); + dim20_ = (mwSize) mxWrapGetScalar(prhs[20], &mw_err_txt_); + dim21_ = (mwSize) mxWrapGetScalar(prhs[21], &mw_err_txt_); + dim22_ = (mwSize) mxWrapGetScalar(prhs[22], &mw_err_txt_); + dim23_ = (mwSize) mxWrapGetScalar(prhs[23], &mw_err_txt_); + dim24_ = (mwSize) mxWrapGetScalar(prhs[24], &mw_err_txt_); + dim25_ = (mwSize) mxWrapGetScalar(prhs[25], &mw_err_txt_); + + if (mxGetM(prhs[0])*mxGetN(prhs[0]) != dim10_) { + mw_err_txt_ = "Bad argument size: nd"; goto mw_err_label; + } + + if (mxGetM(prhs[1]) != dim11_ || + mxGetN(prhs[1]) != dim12_) { + mw_err_txt_ = "Bad argument size: sources"; + goto mw_err_label; + } + + if (mxGetM(prhs[2]) != dim13_ || + mxGetN(prhs[2]) != dim14_) { + mw_err_txt_ = "Bad argument size: stoklet"; + goto mw_err_label; + } + + if (mxGetM(prhs[3])*mxGetN(prhs[3]) != dim15_) { + mw_err_txt_ = "Bad argument size: ns"; goto mw_err_label; + } + + if (mxGetM(prhs[4]) != dim16_ || + mxGetN(prhs[4]) != dim17_) { + mw_err_txt_ = "Bad argument size: targ"; + goto mw_err_label; + } + + if (mxGetM(prhs[5])*mxGetN(prhs[5]) != dim18_) { + mw_err_txt_ = "Bad argument size: nt"; goto mw_err_label; + } + + if (mxGetM(prhs[6]) != dim19_ || + mxGetN(prhs[6]) != dim20_) { + mw_err_txt_ = "Bad argument size: pottarg"; + goto mw_err_label; + } + + if (mxGetM(prhs[7]) != dim21_ || + mxGetN(prhs[7]) != dim22_) { + mw_err_txt_ = "Bad argument size: pretarg"; + goto mw_err_label; + } + + if (mxGetM(prhs[8]) != dim23_ || + mxGetN(prhs[8]) != dim24_) { + mw_err_txt_ = "Bad argument size: gradtarg"; + goto mw_err_label; + } + + if (mxGetM(prhs[9])*mxGetN(prhs[9]) != dim25_) { + mw_err_txt_ = "Bad argument size: thresh"; goto mw_err_label; + } + + if (mxGetM(prhs[0])*mxGetN(prhs[0]) != 0) { + in0_ = mxWrapGetArray_int(prhs[0], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in0_ = NULL; + if (mxGetM(prhs[1])*mxGetN(prhs[1]) != 0) { + if( mxGetClassID(prhs[1]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in1_ = mxGetDoubles(prhs[1]); +#else + in1_ = mxGetPr(prhs[1]); +#endif + } else + in1_ = NULL; + if (mxGetM(prhs[2])*mxGetN(prhs[2]) != 0) { + if( mxGetClassID(prhs[2]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in2_ = mxGetDoubles(prhs[2]); +#else + in2_ = mxGetPr(prhs[2]); +#endif + } else + in2_ = NULL; + if (mxGetM(prhs[3])*mxGetN(prhs[3]) != 0) { + in3_ = mxWrapGetArray_int(prhs[3], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in3_ = NULL; + if (mxGetM(prhs[4])*mxGetN(prhs[4]) != 0) { + if( mxGetClassID(prhs[4]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in4_ = mxGetDoubles(prhs[4]); +#else + in4_ = mxGetPr(prhs[4]); +#endif + } else + in4_ = NULL; + if (mxGetM(prhs[5])*mxGetN(prhs[5]) != 0) { + in5_ = mxWrapGetArray_int(prhs[5], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in5_ = NULL; + if (mxGetM(prhs[6])*mxGetN(prhs[6]) != 0) { + in6_ = mxWrapGetArray_double(prhs[6], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in6_ = NULL; + if (mxGetM(prhs[7])*mxGetN(prhs[7]) != 0) { + in7_ = mxWrapGetArray_double(prhs[7], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in7_ = NULL; + if (mxGetM(prhs[8])*mxGetN(prhs[8]) != 0) { + in8_ = mxWrapGetArray_double(prhs[8], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in8_ = NULL; + if (mxGetM(prhs[9])*mxGetN(prhs[9]) != 0) { + if( mxGetClassID(prhs[9]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in9_ = mxGetDoubles(prhs[9]); +#else + in9_ = mxGetPr(prhs[9]); +#endif + } else + in9_ = NULL; + if (mexprofrecord_) + mexprofrecord_[46]++; + MWF77_st2ddirectstokg(in0_, in1_, in2_, in3_, in4_, in5_, in6_, in7_, in8_, in9_); + plhs[0] = mxCreateDoubleMatrix(dim19_, dim20_, mxREAL); + mxWrapCopy_double(plhs[0], in6_, dim19_*dim20_); + plhs[1] = mxCreateDoubleMatrix(dim21_, dim22_, mxREAL); + mxWrapCopy_double(plhs[1], in7_, dim21_*dim22_); + plhs[2] = mxCreateDoubleMatrix(dim23_, dim24_, mxREAL); + mxWrapCopy_double(plhs[2], in8_, dim23_*dim24_); + +mw_err_label: + if (in0_) mxFree(in0_); + if (in3_) mxFree(in3_); + if (in5_) mxFree(in5_); + if (in6_) mxFree(in6_); + if (in7_) mxFree(in7_); + if (in8_) mxFree(in8_); + if (mw_err_txt_) + mexErrMsgTxt(mw_err_txt_); +} + +/* ---- fmm2d.mw: 1747 ---- + * st2ddirectstokstrsg(int[1] nd, double[2, ns] sources, int[1] ifstoklet, double[nd2, ns_stok] stoklet, int[1] istress, double[nd2, ns_strs] strslet, double[nd2, ns_strs] strsvec, int[1] ns, double[2, nt] targ, int[1] nt, inout double[nd2, nt] pottarg, inout double[nd, nt] pretarg, inout double[nd4, nt] gradtarg, double[1] thresh); + */ +static const char* stubids47_ = "st2ddirectstokstrsg(i int[x], i double[xx], i int[x], i double[xx], i int[x], i double[xx], i double[xx], i int[x], i double[xx], i int[x], io double[xx], io double[xx], io double[xx], i double[x])"; + +void mexStub47(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + const char* mw_err_txt_ = 0; + int* in0_ =0; /* nd */ + double* in1_ =0; /* sources */ + int* in2_ =0; /* ifstoklet */ + double* in3_ =0; /* stoklet */ + int* in4_ =0; /* istress */ + double* in5_ =0; /* strslet */ + double* in6_ =0; /* strsvec */ + int* in7_ =0; /* ns */ + double* in8_ =0; /* targ */ + int* in9_ =0; /* nt */ + double* in10_ =0; /* pottarg */ + double* in11_ =0; /* pretarg */ + double* in12_ =0; /* gradtarg */ + double* in13_ =0; /* thresh */ + mwSize dim14_; /* 1 */ + mwSize dim15_; /* 2 */ + mwSize dim16_; /* ns */ + mwSize dim17_; /* 1 */ + mwSize dim18_; /* nd2 */ + mwSize dim19_; /* ns_stok */ + mwSize dim20_; /* 1 */ + mwSize dim21_; /* nd2 */ + mwSize dim22_; /* ns_strs */ + mwSize dim23_; /* nd2 */ + mwSize dim24_; /* ns_strs */ + mwSize dim25_; /* 1 */ + mwSize dim26_; /* 2 */ + mwSize dim27_; /* nt */ + mwSize dim28_; /* 1 */ + mwSize dim29_; /* nd2 */ + mwSize dim30_; /* nt */ + mwSize dim31_; /* nd */ + mwSize dim32_; /* nt */ + mwSize dim33_; /* nd4 */ + mwSize dim34_; /* nt */ + mwSize dim35_; /* 1 */ + + dim14_ = (mwSize) mxWrapGetScalar(prhs[14], &mw_err_txt_); + dim15_ = (mwSize) mxWrapGetScalar(prhs[15], &mw_err_txt_); + dim16_ = (mwSize) mxWrapGetScalar(prhs[16], &mw_err_txt_); + dim17_ = (mwSize) mxWrapGetScalar(prhs[17], &mw_err_txt_); + dim18_ = (mwSize) mxWrapGetScalar(prhs[18], &mw_err_txt_); + dim19_ = (mwSize) mxWrapGetScalar(prhs[19], &mw_err_txt_); + dim20_ = (mwSize) mxWrapGetScalar(prhs[20], &mw_err_txt_); + dim21_ = (mwSize) mxWrapGetScalar(prhs[21], &mw_err_txt_); + dim22_ = (mwSize) mxWrapGetScalar(prhs[22], &mw_err_txt_); + dim23_ = (mwSize) mxWrapGetScalar(prhs[23], &mw_err_txt_); + dim24_ = (mwSize) mxWrapGetScalar(prhs[24], &mw_err_txt_); + dim25_ = (mwSize) mxWrapGetScalar(prhs[25], &mw_err_txt_); + dim26_ = (mwSize) mxWrapGetScalar(prhs[26], &mw_err_txt_); + dim27_ = (mwSize) mxWrapGetScalar(prhs[27], &mw_err_txt_); + dim28_ = (mwSize) mxWrapGetScalar(prhs[28], &mw_err_txt_); + dim29_ = (mwSize) mxWrapGetScalar(prhs[29], &mw_err_txt_); + dim30_ = (mwSize) mxWrapGetScalar(prhs[30], &mw_err_txt_); + dim31_ = (mwSize) mxWrapGetScalar(prhs[31], &mw_err_txt_); + dim32_ = (mwSize) mxWrapGetScalar(prhs[32], &mw_err_txt_); + dim33_ = (mwSize) mxWrapGetScalar(prhs[33], &mw_err_txt_); + dim34_ = (mwSize) mxWrapGetScalar(prhs[34], &mw_err_txt_); + dim35_ = (mwSize) mxWrapGetScalar(prhs[35], &mw_err_txt_); + + if (mxGetM(prhs[0])*mxGetN(prhs[0]) != dim14_) { + mw_err_txt_ = "Bad argument size: nd"; goto mw_err_label; + } + + if (mxGetM(prhs[1]) != dim15_ || + mxGetN(prhs[1]) != dim16_) { + mw_err_txt_ = "Bad argument size: sources"; + goto mw_err_label; + } + + if (mxGetM(prhs[2])*mxGetN(prhs[2]) != dim17_) { + mw_err_txt_ = "Bad argument size: ifstoklet"; goto mw_err_label; + } + + if (mxGetM(prhs[3]) != dim18_ || + mxGetN(prhs[3]) != dim19_) { + mw_err_txt_ = "Bad argument size: stoklet"; + goto mw_err_label; + } + + if (mxGetM(prhs[4])*mxGetN(prhs[4]) != dim20_) { + mw_err_txt_ = "Bad argument size: istress"; goto mw_err_label; + } + + if (mxGetM(prhs[5]) != dim21_ || + mxGetN(prhs[5]) != dim22_) { + mw_err_txt_ = "Bad argument size: strslet"; + goto mw_err_label; + } + + if (mxGetM(prhs[6]) != dim23_ || + mxGetN(prhs[6]) != dim24_) { + mw_err_txt_ = "Bad argument size: strsvec"; + goto mw_err_label; + } + + if (mxGetM(prhs[7])*mxGetN(prhs[7]) != dim25_) { + mw_err_txt_ = "Bad argument size: ns"; goto mw_err_label; + } + + if (mxGetM(prhs[8]) != dim26_ || + mxGetN(prhs[8]) != dim27_) { + mw_err_txt_ = "Bad argument size: targ"; + goto mw_err_label; + } + + if (mxGetM(prhs[9])*mxGetN(prhs[9]) != dim28_) { + mw_err_txt_ = "Bad argument size: nt"; goto mw_err_label; + } + + if (mxGetM(prhs[10]) != dim29_ || + mxGetN(prhs[10]) != dim30_) { + mw_err_txt_ = "Bad argument size: pottarg"; + goto mw_err_label; + } + + if (mxGetM(prhs[11]) != dim31_ || + mxGetN(prhs[11]) != dim32_) { + mw_err_txt_ = "Bad argument size: pretarg"; + goto mw_err_label; + } + + if (mxGetM(prhs[12]) != dim33_ || + mxGetN(prhs[12]) != dim34_) { + mw_err_txt_ = "Bad argument size: gradtarg"; + goto mw_err_label; + } + + if (mxGetM(prhs[13])*mxGetN(prhs[13]) != dim35_) { + mw_err_txt_ = "Bad argument size: thresh"; goto mw_err_label; + } + + if (mxGetM(prhs[0])*mxGetN(prhs[0]) != 0) { + in0_ = mxWrapGetArray_int(prhs[0], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in0_ = NULL; + if (mxGetM(prhs[1])*mxGetN(prhs[1]) != 0) { + if( mxGetClassID(prhs[1]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in1_ = mxGetDoubles(prhs[1]); +#else + in1_ = mxGetPr(prhs[1]); +#endif + } else + in1_ = NULL; + if (mxGetM(prhs[2])*mxGetN(prhs[2]) != 0) { + in2_ = mxWrapGetArray_int(prhs[2], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in2_ = NULL; + if (mxGetM(prhs[3])*mxGetN(prhs[3]) != 0) { + if( mxGetClassID(prhs[3]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in3_ = mxGetDoubles(prhs[3]); +#else + in3_ = mxGetPr(prhs[3]); +#endif + } else + in3_ = NULL; + if (mxGetM(prhs[4])*mxGetN(prhs[4]) != 0) { + in4_ = mxWrapGetArray_int(prhs[4], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in4_ = NULL; + if (mxGetM(prhs[5])*mxGetN(prhs[5]) != 0) { + if( mxGetClassID(prhs[5]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in5_ = mxGetDoubles(prhs[5]); +#else + in5_ = mxGetPr(prhs[5]); +#endif + } else + in5_ = NULL; + if (mxGetM(prhs[6])*mxGetN(prhs[6]) != 0) { + if( mxGetClassID(prhs[6]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in6_ = mxGetDoubles(prhs[6]); +#else + in6_ = mxGetPr(prhs[6]); +#endif + } else + in6_ = NULL; + if (mxGetM(prhs[7])*mxGetN(prhs[7]) != 0) { + in7_ = mxWrapGetArray_int(prhs[7], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in7_ = NULL; + if (mxGetM(prhs[8])*mxGetN(prhs[8]) != 0) { + if( mxGetClassID(prhs[8]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in8_ = mxGetDoubles(prhs[8]); +#else + in8_ = mxGetPr(prhs[8]); +#endif + } else + in8_ = NULL; + if (mxGetM(prhs[9])*mxGetN(prhs[9]) != 0) { + in9_ = mxWrapGetArray_int(prhs[9], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in9_ = NULL; + if (mxGetM(prhs[10])*mxGetN(prhs[10]) != 0) { + in10_ = mxWrapGetArray_double(prhs[10], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in10_ = NULL; + if (mxGetM(prhs[11])*mxGetN(prhs[11]) != 0) { + in11_ = mxWrapGetArray_double(prhs[11], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in11_ = NULL; + if (mxGetM(prhs[12])*mxGetN(prhs[12]) != 0) { + in12_ = mxWrapGetArray_double(prhs[12], &mw_err_txt_); + if (mw_err_txt_) + goto mw_err_label; + } else + in12_ = NULL; + if (mxGetM(prhs[13])*mxGetN(prhs[13]) != 0) { + if( mxGetClassID(prhs[13]) != mxDOUBLE_CLASS ) + mw_err_txt_ = "Invalid array argument, mxDOUBLE_CLASS expected"; + if (mw_err_txt_) goto mw_err_label; +#if MX_HAS_INTERLEAVED_COMPLEX + in13_ = mxGetDoubles(prhs[13]); +#else + in13_ = mxGetPr(prhs[13]); +#endif + } else + in13_ = NULL; + if (mexprofrecord_) + mexprofrecord_[47]++; + MWF77_st2ddirectstokstrsg(in0_, in1_, in2_, in3_, in4_, in5_, in6_, in7_, in8_, in9_, in10_, in11_, in12_, in13_); + plhs[0] = mxCreateDoubleMatrix(dim29_, dim30_, mxREAL); + mxWrapCopy_double(plhs[0], in10_, dim29_*dim30_); + plhs[1] = mxCreateDoubleMatrix(dim31_, dim32_, mxREAL); + mxWrapCopy_double(plhs[1], in11_, dim31_*dim32_); + plhs[2] = mxCreateDoubleMatrix(dim33_, dim34_, mxREAL); + mxWrapCopy_double(plhs[2], in12_, dim33_*dim34_); + +mw_err_label: + if (in0_) mxFree(in0_); + if (in2_) mxFree(in2_); + if (in4_) mxFree(in4_); + if (in7_) mxFree(in7_); + if (in9_) mxFree(in9_); + if (in10_) mxFree(in10_); + if (in11_) mxFree(in11_); + if (in12_) mxFree(in12_); + if (mw_err_txt_) + mexErrMsgTxt(mw_err_txt_); +} + /* ---- */ void mexFunction(int nlhs, mxArray* plhs[], @@ -11186,12 +12067,18 @@ void mexFunction(int nlhs, mxArray* plhs[], mexStub43(nlhs,plhs, nrhs-1,prhs+1); else if (strcmp(id, stubids44_) == 0) mexStub44(nlhs,plhs, nrhs-1,prhs+1); + else if (strcmp(id, stubids45_) == 0) + mexStub45(nlhs,plhs, nrhs-1,prhs+1); + else if (strcmp(id, stubids46_) == 0) + mexStub46(nlhs,plhs, nrhs-1,prhs+1); + else if (strcmp(id, stubids47_) == 0) + mexStub47(nlhs,plhs, nrhs-1,prhs+1); else if (strcmp(id, "*profile on*") == 0) { if (!mexprofrecord_) { - mexprofrecord_ = (int*) malloc(45 * sizeof(int)); + mexprofrecord_ = (int*) malloc(48 * sizeof(int)); mexLock(); } - memset(mexprofrecord_, 0, 45 * sizeof(int)); + memset(mexprofrecord_, 0, 48 * sizeof(int)); } else if (strcmp(id, "*profile off*") == 0) { if (mexprofrecord_) { free(mexprofrecord_); @@ -11243,6 +12130,9 @@ void mexFunction(int nlhs, mxArray* plhs[], mexPrintf("%d calls to fmm2d.mw:1451\n", mexprofrecord_[42]); mexPrintf("%d calls to fmm2d.mw:1454\n", mexprofrecord_[43]); mexPrintf("%d calls to fmm2d.mw:1457\n", mexprofrecord_[44]); + mexPrintf("%d calls to fmm2d.mw:1647\n", mexprofrecord_[45]); + mexPrintf("%d calls to fmm2d.mw:1744\n", mexprofrecord_[46]); + mexPrintf("%d calls to fmm2d.mw:1747\n", mexprofrecord_[47]); } else if (strcmp(id, "*profile log*") == 0) { FILE* logfp; if (nrhs != 2 || mxGetString(prhs[1], id, sizeof(id)) != 0) @@ -11294,6 +12184,9 @@ void mexFunction(int nlhs, mxArray* plhs[], fprintf(logfp, "%d calls to fmm2d.mw:1451\n", mexprofrecord_[42]); fprintf(logfp, "%d calls to fmm2d.mw:1454\n", mexprofrecord_[43]); fprintf(logfp, "%d calls to fmm2d.mw:1457\n", mexprofrecord_[44]); + fprintf(logfp, "%d calls to fmm2d.mw:1647\n", mexprofrecord_[45]); + fprintf(logfp, "%d calls to fmm2d.mw:1744\n", mexprofrecord_[46]); + fprintf(logfp, "%d calls to fmm2d.mw:1747\n", mexprofrecord_[47]); fclose(logfp); } else mexErrMsgTxt("Unknown identifier"); diff --git a/matlab/fmm2d.mw b/matlab/fmm2d.mw index 0697cfe..1d655ec 100644 --- a/matlab/fmm2d.mw +++ b/matlab/fmm2d.mw @@ -1464,3 +1464,293 @@ end % % % --------------------------------------------------------------------- +@function [U] = stfmm2d(eps,srcinfo,ifppreg,targ,ifppregtarg) +% STFMM2d FMM in 2d for Stokes (viscous fluid hydrodynamic) kernels. +% +% U = stfmm2d(eps,srcinfo,ifppreg,targ,ifppregtarg) +% +% Stokes FMM in R^2: evaluate all pairwise particle +% interactions (ignoring self-interactions) and, possibly, +% interactions with targets, to precision eps, using the fast +% multipole method (linear scaling in number of points). +% +% This routine computes the sum for the velocity vector, +% +% u_i(x) = sum_m G_{ij}(x,y^{(m)}) sigma^{(m)}_j +% + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k +% +% for each requested evaluation point x and i=1,2, +% where sigma^{(m)} is the Stokeslet force, mu^{(m)} is the +% stresslet strength, and nu^{(m)} is the stresslet orientation +% (note that each of these is a 2 vector per source point y^{(m)}). +% The stokeslet kernel G and stresslet kernel T are defined below. +% Repeated indices are taken as summed over 1,2, ie, Einstein +% convention. +% +% There are two options for the evaluation points $x$. These +% can be the sources themselves (ifppreg > 0; see below) or other +% target points of interest (ifppregtarg > 0; see below). Both options +% can be selected in one call. +% When $x=y_{(m)}$, the term corresponding to $y^{(m)}$ is dropped +% from the sum. +% +% Optionally, the associated pressure p(x) and 2x2 gradient tensor +% grad u(x) are returned, +% +% p(x) = sum_m P_j(x,y^m) sigma^{(m)}_j +% + sum_m PI_{jk}(x,y{(m)}) mu^{(m)}_j nu^{(m)}_k +% +% grad_l u_i(x) = grad_l [sum_m G_{ij}(x,y^m) sigma^{(m)}_j +% + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k] +% +% where the pressure stokeslet P and pressure stresslet PI are defined +% below. Note that these two may be combined to get the stress tensor. +% +% We use the following as the kernel definitions, noting that: +% 1) The dynamic viscosity (mu) is assumed to be 1. +% 2) All kernels are a factor 2pi larger than standard definitions +% (this is for historical reasons). +% Thus, in general, divide the velocity (potential or grad) outputs +% by i2pi.mu, and pressure by 2pi, to recover standard definitions. +% +% For a source y and target x, let r_i = x_i-y_i (note sign) +% and let r = sqrt(r_1^2 + r_2^2) +% +% The Stokeslet, G_{ij}, and its associated pressure tensor, P_j, +% we define as +% +% G_{ij}(x,y) = (-delta_{ij}log(r) + r_i r_j / r^2 )/2 +% P_j(x,y) = r_j/r^2 +% +% The (Type I) stresslet, T_{ijk}, and its associated pressure +% tensor, PI_{jk}, we define as +% +% T_{ijk}(x,y) = -2 r_i r_j r_k/ r^4 +% PI_{jk}(x,y) = -2 delta_{jk}/r^2 + 2 r_j r_k/r^4 +% +% Args: +% +% - eps: double +% precision requested +% - srcinfo: structure +% structure containing the following info about sources: +% * srcinfo.sources: double(2,n) +% source locations, $x_{j}$ +% * srcinfo.nd: integer +% number of density vectors for the same points (optional, +% default - nd = 1) +% * srcinfo.stoklet: double(nd,2,n) +% Stokeslet charge strengths, $sigma_{j}$ (optional, +% default - term corresponding to Stokeslet charge strengths dropped) +% * srcinfo.strslet: double(nd,2,n) +% stresslet strengths, $mu_{j}$ (optional +% default - term corresponding to stresslet strengths dropped) +% * srcinfo.strsvec: double(nd,2,n) +% stresslet orientations, $nu_{j}$ (optional +% default - term corresponding to stresslet orientations dropped) +% - ifppreg: integer +% | source eval flag +% | potential at sources evaluated if ifppreg = 1 +% | potential and pressure at sources evaluated if ifppreg=2 +% | potential, pressure and gradient at sources evaluated if ifppreg=3 +% - targ: double(2,nt) +% target locations, $t_{i}$ (optional) +% - ifppregtarg: integer +% | target eval flag (optional) +% | potential at targets evaluated if ifppregtarg = 1 +% | potential and pressure at targets evaluated if ifppregtarg = 2 +% | potential, pressure and gradient at targets evaluated if +% | ifppregtarg = 3 +% +% Returns: +% +% - U.pot: velocity at source locations if requested +% - U.pre: pressure at source locations if requested +% - U.grad: gradient of velocity at source locations if requested +% - U.pottarg: velocity at target locations if requested +% - U.pretarg: pressure at target locations if requested +% - U.gradtarg: gradient of velocity at target locations if requested +% +% See also: ST2dDIR + + sources = srcinfo.sources; + [m,ns] = size(sources); + assert(m==2,'The first dimension of sources must be 2'); + if(~isfield(srcinfo,'nd')) + nd = 1; + end + if(isfield(srcinfo,'nd')) + nd = srcinfo.nd; + end + + pot = zeros(nd*2,1); ns_pot = 1; + pre = zeros(nd,1); ns_pre = 1; + grad = zeros(nd*4,1); ns_grad = 1; + + if(ifppreg >= 1), pot = zeros(nd*2,ns); ns_pot = ns; end; + if(ifppreg >= 2), pre = zeros(nd,ns); ns_pre = ns; end; + if(ifppreg >= 3), grad = zeros(nd*4,ns); ns_grad = ns; end; + + pottarg = zeros(nd*2,1); nt_pot = 1; + pretarg = zeros(nd,1); nt_pre = 1; + gradtarg = zeros(nd*4,1); nt_grad = 1; + if( nargin <= 3 ) + nt = 0; + ifppregtarg = 0; + targ = zeros(2,0); + else + if( nargin <= 4 ), ifppregtarg = 0; end; + [m,nt] = size(targ); + assert(m==2,'First dimension of targets must be 2'); + if(ifppregtarg >= 1), pottarg = zeros(nd*2,nt); nt_pot = nt; end; + if(ifppregtarg >= 2), pretarg = zeros(nd,nt); nt_pre = nt; end; + if(ifppregtarg >= 3), gradtarg = zeros(nd*4,nt); nt_grad = nt; end; + end + + if(ifppreg ==0 && ifppregtarg ==0), disp('Nothing to compute, set eigher ifppreg or ifppregtarg to 1 or 2 or 3'); return; end; + + if(isfield(srcinfo,'stoklet')) + ifstoklet = 1; + ns_stok = ns; + stoklet = srcinfo.stoklet; + if(nd == 1), [a,b] = size(squeeze(stoklet)); assert(a==2 && b==ns,'Stoklet must be of shape[2,ns], where ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(stoklet); assert(a==nd && b==2 && c==ns, 'Stoklet must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + stoklet = reshape(stoklet,[2*nd,ns]); + else + ifstoklet = 0; + ns_stok = 1; + stoklet = zeros(nd*2,1); + end + + if(isfield(srcinfo,'strslet') && isfield(srcinfo,'strsvec')) + ifstrslet = 1; + ns_strs = ns; + strslet = srcinfo.strslet; + strsvec = srcinfo.strsvec; + if(nd == 1), [a,b] = size(squeeze(strslet)); assert(a==2 && b==ns,'Strslet must be of shape[2,ns], where ns is the number of sources'); end; + if(nd == 1), [a,b] = size(squeeze(strsvec)); assert(a==2 && b==ns,'Strsvec must be of shape[2,ns], where ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(strslet); assert(a==nd && b==2 && c==ns, 'Strslet must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(strsvec); assert(a==nd && b==2 && c==ns, 'Strsvec must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + strslet = reshape(strslet,[2*nd,ns]); + strsvec = reshape(strsvec,[2*nd,ns]); + else + ifstrslet = 0; + ns_strs = 1; + strslet = zeros(nd*2,1); + strsvec = zeros(nd*2,1); + end + + nd2 = 2*nd; + nd4 = 4*nd; + ier = 0; + + # FORTRAN stfmm2d(int[1] nd, double[1] eps, int[1] ns, double[2,ns] sources, int[1] ifstoklet, double[nd2,ns_stok] stoklet, int[1] ifstrslet, double[nd2,ns_strs] strslet, double[nd2,ns_strs] strsvec, int[1] ifppreg, inout double[nd2,ns_pot] pot, inout double[nd,ns_pre] pre, inout double[nd4,ns_grad] grad, int[1] nt, double[2,nt] targ, int[1] ifppregtarg, inout double[nd2,nt_pot] pottarg, inout double [nd,nt_pre] pretarg, inout double[nd4,nt_grad] gradtarg, inout int[1] ier); + + U.pot = []; + U.pre = []; + U.grad = []; + U.pottarg = []; + U.pretarg = []; + U.gradtarg = []; + if(ifppreg >= 1), U.pot = squeeze(reshape(pot,[nd,2,ns])); end; + if(ifppreg >= 2), U.pre = pre; end; + if(ifppreg >= 3), U.grad = squeeze(reshape(grad,[nd,2,2,ns])); end; + if(ifppregtarg >= 1), U.pottarg = squeeze(reshape(pottarg,[nd,2,nt])); end; + if(ifppregtarg >= 2), U.pretarg = pretarg; end; + if(ifppregtarg >= 3), U.gradtarg = squeeze(reshape(gradtarg,[nd,2,2,nt])); end; + +end + +% --------------------------------------------------------------------- +@function [U] = st2ddir(srcinfo,targ,ifppregtarg) +% ST2dDIR Direct (slow) 2d Stokes kernel sums (reference for STFMM2d). +% +% U = st2ddir(srcinfo,targ,ifppregtarg) +% +% Stokes direct evaluation in R^2: evaluate all pairwise particle +% interactions with targets. This is the slow O(N^2) direct code used +% as a reference for testing the (fast) code stfmm2d. +% +% Kernel definitions, input and outputs arguments are identical to +% stfmm2d (see that function for all definitions), apart from: +% 1) the first argument (eps) is absent. +% 2) there are currently no outputs at sources, meaning that U.pot, U.pre, +% and U.grad are missing (as if ifppreg=0). In other words, +% just targets for now, and targ is thus not an optional argument. +% +% See also: STFMM2d + + sources = srcinfo.sources; + [m,ns] = size(sources); + assert(m==2,'The first dimension of sources must be 2'); + if(~isfield(srcinfo,'nd')) + nd = 1; + end + if(isfield(srcinfo,'nd')) + nd = srcinfo.nd; + end + + thresh = 1e-15; + + if( nargin <= 1 ) + return; + else + if( nargin <= 2 ), ifppregtarg = 3; end; + [m,nt] = size(targ); + assert(m==2,'First dimension of targets must be 2'); + pottarg = zeros(nd*2,nt); + pretarg = zeros(nd,nt); + gradtarg = zeros(nd*4,nt); + end + + if(ifppregtarg == 0), disp('Nothing to compute, set eigher ifppregtarg to 1 or 2 or 3'); return; end; + + if(isfield(srcinfo,'stoklet')) + ifstoklet = 1; + ns_stok = ns; + stoklet = srcinfo.stoklet; + if(nd == 1), [a,b] = size(squeeze(stoklet)); assert(a==2 && b==ns,'Stoklet must be of shape[2,ns], where ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(stoklet); assert(a==nd && b==2 && c==ns, 'Stoklet must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + stoklet = reshape(stoklet,[2*nd,ns]); + else + ifstoklet = 0; + ns_stok = 1; + stoklet = zeros(nd*2,1); + end + + if(isfield(srcinfo,'strslet') && isfield(srcinfo,'strsvec')) + ifstrslet = 1; + ns_strs = ns; + strslet = srcinfo.strslet; + strsvec = srcinfo.strsvec; + if(nd == 1), [a,b] = size(squeeze(strslet)); assert(a==2 && b==ns,'Strslet must be of shape[2,ns], where ns is the number of sources'); end; + if(nd == 1), [a,b] = size(squeeze(strsvec)); assert(a==2 && b==ns,'Strsvec must be of shape[2,ns], where ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(strslet); assert(a==nd && b==2 && c==ns, 'Strslet must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(strsvec); assert(a==nd && b==2 && c==ns, 'Strsvec must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + strslet = reshape(strslet,[2*nd,ns]); + strsvec = reshape(strsvec,[2*nd,ns]); + else + ifstrslet = 0; + ns_strs = 1; + strslet = zeros(nd*2,1); + strsvec = zeros(nd*2,1); + end + + nd2 = 2*nd; + nd4 = 4*nd; + ier = 0; + + if ifstoklet == 1 && ifstrslet == 0 + # FORTRAN st2ddirectstokg(int[1] nd, double[2,ns] sources, double[nd2,ns_stok] stoklet, int[1] ns, double[2,nt] targ, int[1] nt, inout double[nd2,nt] pottarg, inout double [nd,nt] pretarg, inout double[nd4,nt] gradtarg, double[1] thresh); + else + istress = 1; + # FORTRAN st2ddirectstokstrsg(int[1] nd, double[2,ns] sources, int[1] ifstoklet, double[nd2,ns_stok] stoklet, int[1] istress, double[nd2,ns_strs] strslet, double[nd2,ns_strs] strsvec, int[1] ns, double[2,nt] targ, int[1] nt, inout double[nd2,nt] pottarg, inout double [nd,nt] pretarg, inout double[nd4,nt] gradtarg, double[1] thresh); + end + + U.pottarg = []; + U.pretarg = []; + U.gradtarg = []; + if(ifppregtarg >= 1), U.pottarg = squeeze(reshape(pottarg,[nd,2,nt])); end; + if(ifppregtarg >= 2), U.pretarg = pretarg; end; + if(ifppregtarg >= 3), U.gradtarg = squeeze(reshape(gradtarg,[nd,2,2,nt])); end; +end diff --git a/matlab/st2ddir.m b/matlab/st2ddir.m new file mode 100644 index 0000000..029a293 --- /dev/null +++ b/matlab/st2ddir.m @@ -0,0 +1,94 @@ +function [U] = st2ddir(srcinfo,targ,ifppregtarg) +% ST2dDIR Direct (slow) 2d Stokes kernel sums (reference for STFMM2d). +% +% U = st2ddir(srcinfo,targ,ifppregtarg) +% +% Stokes direct evaluation in R^2: evaluate all pairwise particle +% interactions with targets. This is the slow O(N^2) direct code used +% as a reference for testing the (fast) code stfmm2d. +% +% Kernel definitions, input and outputs arguments are identical to +% stfmm2d (see that function for all definitions), apart from: +% 1) the first argument (eps) is absent. +% 2) there are currently no outputs at sources, meaning that U.pot, U.pre, +% and U.grad are missing (as if ifppreg=0). In other words, +% just targets for now, and targ is thus not an optional argument. +% +% See also: STFMM2d + + sources = srcinfo.sources; + [m,ns] = size(sources); + assert(m==2,'The first dimension of sources must be 2'); + if(~isfield(srcinfo,'nd')) + nd = 1; + end + if(isfield(srcinfo,'nd')) + nd = srcinfo.nd; + end + + thresh = 1e-15; + + if( nargin <= 1 ) + return; + else + if( nargin <= 2 ), ifppregtarg = 3; end; + [m,nt] = size(targ); + assert(m==2,'First dimension of targets must be 2'); + pottarg = zeros(nd*2,nt); + pretarg = zeros(nd,nt); + gradtarg = zeros(nd*4,nt); + end + + if(ifppregtarg == 0), disp('Nothing to compute, set eigher ifppregtarg to 1 or 2 or 3'); return; end; + + if(isfield(srcinfo,'stoklet')) + ifstoklet = 1; + ns_stok = ns; + stoklet = srcinfo.stoklet; + if(nd == 1), [a,b] = size(squeeze(stoklet)); assert(a==2 && b==ns,'Stoklet must be of shape[2,ns], where ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(stoklet); assert(a==nd && b==2 && c==ns, 'Stoklet must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + stoklet = reshape(stoklet,[2*nd,ns]); + else + ifstoklet = 0; + ns_stok = 1; + stoklet = zeros(nd*2,1); + end + + if(isfield(srcinfo,'strslet') && isfield(srcinfo,'strsvec')) + ifstrslet = 1; + ns_strs = ns; + strslet = srcinfo.strslet; + strsvec = srcinfo.strsvec; + if(nd == 1), [a,b] = size(squeeze(strslet)); assert(a==2 && b==ns,'Strslet must be of shape[2,ns], where ns is the number of sources'); end; + if(nd == 1), [a,b] = size(squeeze(strsvec)); assert(a==2 && b==ns,'Strsvec must be of shape[2,ns], where ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(strslet); assert(a==nd && b==2 && c==ns, 'Strslet must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(strsvec); assert(a==nd && b==2 && c==ns, 'Strsvec must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + strslet = reshape(strslet,[2*nd,ns]); + strsvec = reshape(strsvec,[2*nd,ns]); + else + ifstrslet = 0; + ns_strs = 1; + strslet = zeros(nd*2,1); + strsvec = zeros(nd*2,1); + end + + nd2 = 2*nd; + nd4 = 4*nd; + ier = 0; + + if ifstoklet == 1 && ifstrslet == 0 + mex_id_ = 'st2ddirectstokg(i int[x], i double[xx], i double[xx], i int[x], i double[xx], i int[x], io double[xx], io double[xx], io double[xx], i double[x])'; +[pottarg, pretarg, gradtarg] = fmm2d(mex_id_, nd, sources, stoklet, ns, targ, nt, pottarg, pretarg, gradtarg, thresh, 1, 2, ns, nd2, ns_stok, 1, 2, nt, 1, nd2, nt, nd, nt, nd4, nt, 1); + else + istress = 1; + mex_id_ = 'st2ddirectstokstrsg(i int[x], i double[xx], i int[x], i double[xx], i int[x], i double[xx], i double[xx], i int[x], i double[xx], i int[x], io double[xx], io double[xx], io double[xx], i double[x])'; +[pottarg, pretarg, gradtarg] = fmm2d(mex_id_, nd, sources, ifstoklet, stoklet, istress, strslet, strsvec, ns, targ, nt, pottarg, pretarg, gradtarg, thresh, 1, 2, ns, 1, nd2, ns_stok, 1, nd2, ns_strs, nd2, ns_strs, 1, 2, nt, 1, nd2, nt, nd, nt, nd4, nt, 1); + end + + U.pottarg = []; + U.pretarg = []; + U.gradtarg = []; + if(ifppregtarg >= 1), U.pottarg = squeeze(reshape(pottarg,[nd,2,nt])); end; + if(ifppregtarg >= 2), U.pretarg = pretarg; end; + if(ifppregtarg >= 3), U.gradtarg = squeeze(reshape(gradtarg,[nd,2,2,nt])); end; +end diff --git a/matlab/stfmm2d.m b/matlab/stfmm2d.m new file mode 100644 index 0000000..6df3f62 --- /dev/null +++ b/matlab/stfmm2d.m @@ -0,0 +1,199 @@ +function [U] = stfmm2d(eps,srcinfo,ifppreg,targ,ifppregtarg) +% STFMM2d FMM in 2d for Stokes (viscous fluid hydrodynamic) kernels. +% +% U = stfmm2d(eps,srcinfo,ifppreg,targ,ifppregtarg) +% +% Stokes FMM in R^2: evaluate all pairwise particle +% interactions (ignoring self-interactions) and, possibly, +% interactions with targets, to precision eps, using the fast +% multipole method (linear scaling in number of points). +% +% This routine computes the sum for the velocity vector, +% +% u_i(x) = sum_m G_{ij}(x,y^{(m)}) sigma^{(m)}_j +% + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k +% +% for each requested evaluation point x and i=1,2, +% where sigma^{(m)} is the Stokeslet force, mu^{(m)} is the +% stresslet strength, and nu^{(m)} is the stresslet orientation +% (note that each of these is a 2 vector per source point y^{(m)}). +% The stokeslet kernel G and stresslet kernel T are defined below. +% Repeated indices are taken as summed over 1,2, ie, Einstein +% convention. +% +% There are two options for the evaluation points $x$. These +% can be the sources themselves (ifppreg > 0; see below) or other +% target points of interest (ifppregtarg > 0; see below). Both options +% can be selected in one call. +% When $x=y_{(m)}$, the term corresponding to $y^{(m)}$ is dropped +% from the sum. +% +% Optionally, the associated pressure p(x) and 2x2 gradient tensor +% grad u(x) are returned, +% +% p(x) = sum_m P_j(x,y^m) sigma^{(m)}_j +% + sum_m PI_{jk}(x,y{(m)}) mu^{(m)}_j nu^{(m)}_k +% +% grad_l u_i(x) = grad_l [sum_m G_{ij}(x,y^m) sigma^{(m)}_j +% + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k] +% +% where the pressure stokeslet P and pressure stresslet PI are defined +% below. Note that these two may be combined to get the stress tensor. +% +% We use the following as the kernel definitions, noting that: +% 1) The dynamic viscosity (mu) is assumed to be 1. +% 2) All kernels are a factor 2pi larger than standard definitions +% (this is for historical reasons). +% Thus, in general, divide the velocity (potential or grad) outputs +% by i2pi.mu, and pressure by 2pi, to recover standard definitions. +% +% For a source y and target x, let r_i = x_i-y_i (note sign) +% and let r = sqrt(r_1^2 + r_2^2) +% +% The Stokeslet, G_{ij}, and its associated pressure tensor, P_j, +% we define as +% +% G_{ij}(x,y) = (-delta_{ij}log(r) + r_i r_j / r^2 )/2 +% P_j(x,y) = r_j/r^2 +% +% The (Type I) stresslet, T_{ijk}, and its associated pressure +% tensor, PI_{jk}, we define as +% +% T_{ijk}(x,y) = -2 r_i r_j r_k/ r^4 +% PI_{jk}(x,y) = -2 delta_{jk}/r^2 + 2 r_j r_k/r^4 +% +% Args: +% +% - eps: double +% precision requested +% - srcinfo: structure +% structure containing the following info about sources: +% * srcinfo.sources: double(2,n) +% source locations, $x_{j}$ +% * srcinfo.nd: integer +% number of density vectors for the same points (optional, +% default - nd = 1) +% * srcinfo.stoklet: double(nd,2,n) +% Stokeslet charge strengths, $sigma_{j}$ (optional, +% default - term corresponding to Stokeslet charge strengths dropped) +% * srcinfo.strslet: double(nd,2,n) +% stresslet strengths, $mu_{j}$ (optional +% default - term corresponding to stresslet strengths dropped) +% * srcinfo.strsvec: double(nd,2,n) +% stresslet orientations, $nu_{j}$ (optional +% default - term corresponding to stresslet orientations dropped) +% - ifppreg: integer +% | source eval flag +% | potential at sources evaluated if ifppreg = 1 +% | potential and pressure at sources evaluated if ifppreg=2 +% | potential, pressure and gradient at sources evaluated if ifppreg=3 +% - targ: double(2,nt) +% target locations, $t_{i}$ (optional) +% - ifppregtarg: integer +% | target eval flag (optional) +% | potential at targets evaluated if ifppregtarg = 1 +% | potential and pressure at targets evaluated if ifppregtarg = 2 +% | potential, pressure and gradient at targets evaluated if +% | ifppregtarg = 3 +% +% Returns: +% +% - U.pot: velocity at source locations if requested +% - U.pre: pressure at source locations if requested +% - U.grad: gradient of velocity at source locations if requested +% - U.pottarg: velocity at target locations if requested +% - U.pretarg: pressure at target locations if requested +% - U.gradtarg: gradient of velocity at target locations if requested +% +% See also: ST2dDIR + + sources = srcinfo.sources; + [m,ns] = size(sources); + assert(m==2,'The first dimension of sources must be 2'); + if(~isfield(srcinfo,'nd')) + nd = 1; + end + if(isfield(srcinfo,'nd')) + nd = srcinfo.nd; + end + + pot = zeros(nd*2,1); ns_pot = 1; + pre = zeros(nd,1); ns_pre = 1; + grad = zeros(nd*4,1); ns_grad = 1; + + if(ifppreg >= 1), pot = zeros(nd*2,ns); ns_pot = ns; end; + if(ifppreg >= 2), pre = zeros(nd,ns); ns_pre = ns; end; + if(ifppreg >= 3), grad = zeros(nd*4,ns); ns_grad = ns; end; + + pottarg = zeros(nd*2,1); nt_pot = 1; + pretarg = zeros(nd,1); nt_pre = 1; + gradtarg = zeros(nd*4,1); nt_grad = 1; + if( nargin <= 3 ) + nt = 0; + ifppregtarg = 0; + targ = zeros(2,0); + else + if( nargin <= 4 ), ifppregtarg = 0; end; + [m,nt] = size(targ); + assert(m==2,'First dimension of targets must be 2'); + if(ifppregtarg >= 1), pottarg = zeros(nd*2,nt); nt_pot = nt; end; + if(ifppregtarg >= 2), pretarg = zeros(nd,nt); nt_pre = nt; end; + if(ifppregtarg >= 3), gradtarg = zeros(nd*4,nt); nt_grad = nt; end; + end + + if(ifppreg ==0 && ifppregtarg ==0), disp('Nothing to compute, set eigher ifppreg or ifppregtarg to 1 or 2 or 3'); return; end; + + if(isfield(srcinfo,'stoklet')) + ifstoklet = 1; + ns_stok = ns; + stoklet = srcinfo.stoklet; + if(nd == 1), [a,b] = size(squeeze(stoklet)); assert(a==2 && b==ns,'Stoklet must be of shape[2,ns], where ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(stoklet); assert(a==nd && b==2 && c==ns, 'Stoklet must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + stoklet = reshape(stoklet,[2*nd,ns]); + else + ifstoklet = 0; + ns_stok = 1; + stoklet = zeros(nd*2,1); + end + + if(isfield(srcinfo,'strslet') && isfield(srcinfo,'strsvec')) + ifstrslet = 1; + ns_strs = ns; + strslet = srcinfo.strslet; + strsvec = srcinfo.strsvec; + if(nd == 1), [a,b] = size(squeeze(strslet)); assert(a==2 && b==ns,'Strslet must be of shape[2,ns], where ns is the number of sources'); end; + if(nd == 1), [a,b] = size(squeeze(strsvec)); assert(a==2 && b==ns,'Strsvec must be of shape[2,ns], where ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(strslet); assert(a==nd && b==2 && c==ns, 'Strslet must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + if(nd>1), [a,b,c] = size(strsvec); assert(a==nd && b==2 && c==ns, 'Strsvec must be of shape[nd,2,ns], where nd is number of densities, and ns is the number of sources'); end; + strslet = reshape(strslet,[2*nd,ns]); + strsvec = reshape(strsvec,[2*nd,ns]); + else + ifstrslet = 0; + ns_strs = 1; + strslet = zeros(nd*2,1); + strsvec = zeros(nd*2,1); + end + + nd2 = 2*nd; + nd4 = 4*nd; + ier = 0; + + mex_id_ = 'stfmm2d(i int[x], i double[x], i int[x], i double[xx], i int[x], i double[xx], i int[x], i double[xx], i double[xx], i int[x], io double[xx], io double[xx], io double[xx], i int[x], i double[xx], i int[x], io double[xx], io double[xx], io double[xx], io int[x])'; +[pot, pre, grad, pottarg, pretarg, gradtarg, ier] = fmm2d(mex_id_, nd, eps, ns, sources, ifstoklet, stoklet, ifstrslet, strslet, strsvec, ifppreg, pot, pre, grad, nt, targ, ifppregtarg, pottarg, pretarg, gradtarg, ier, 1, 1, 1, 2, ns, 1, nd2, ns_stok, 1, nd2, ns_strs, nd2, ns_strs, 1, nd2, ns_pot, nd, ns_pre, nd4, ns_grad, 1, 2, nt, 1, nd2, nt_pot, nd, nt_pre, nd4, nt_grad, 1); + + U.pot = []; + U.pre = []; + U.grad = []; + U.pottarg = []; + U.pretarg = []; + U.gradtarg = []; + if(ifppreg >= 1), U.pot = squeeze(reshape(pot,[nd,2,ns])); end; + if(ifppreg >= 2), U.pre = pre; end; + if(ifppreg >= 3), U.grad = squeeze(reshape(grad,[nd,2,2,ns])); end; + if(ifppregtarg >= 1), U.pottarg = squeeze(reshape(pottarg,[nd,2,nt])); end; + if(ifppregtarg >= 2), U.pretarg = pretarg; end; + if(ifppregtarg >= 3), U.gradtarg = squeeze(reshape(gradtarg,[nd,2,2,nt])); end; + +end + +% --------------------------------------------------------------------- diff --git a/matlab/stfmm2dTest.m b/matlab/stfmm2dTest.m new file mode 100644 index 0000000..3cafc41 --- /dev/null +++ b/matlab/stfmm2dTest.m @@ -0,0 +1,82 @@ +clear srcinfo + +ns = 4000; +srcinfo.sources = rand(2,ns); +srcinfo.stoklet = rand(2,ns); +srcinfo.strslet = rand(2,ns); +srcinfo.strsvec = rand(2,ns); + +nt = 3999; +targ = rand(2,nt); + +eps = 1e-5; +ntests = 2; +ipass = zeros(ntests,1); +errs = zeros(ntests,1); +ntest = 10; + +stmp = srcinfo.sources(:,1:ntest); +ttmp = targ(:,1:ntest); + +itest = 1; + + +% Test source to source, target, stoklet, strslet, pot, pre, grad +ifppreg = 3; +ifppregtarg = 3; +U1 = stfmm2d(eps,srcinfo,ifppreg,targ,ifppregtarg); +U2 = st2ddir(srcinfo,stmp,ifppreg); +U3 = st2ddir(srcinfo,ttmp,ifppregtarg); +s1 = U1.pot(:,1:ntest); s1 = s1(:); +s2 = U2.pottarg; s2 = s2(:); +s3 = U1.pre(1:ntest); s3 = s3(:); +s4 = U2.pretarg; s4 = s4(:); +s5 = U1.grad(:,:,1:ntest); s5 = s5(:); +s6 = U2.gradtarg; s6 = s6(:); +t1 = U1.pottarg(:,1:ntest); t1 = t1(:); +t2 = U3.pottarg; t2 = t2(:); +t3 = U1.pretarg(1:ntest); t3 = t3(:); +t4 = U3.pretarg; t4 = t4(:); +t5 = U1.gradtarg(:,:,1:ntest); t5 = t5(:); +t6 = U3.gradtarg; t6 = t6(:); + +err = norm(s1-s2)^2 + norm(s3-s4)^2 + norm(s5-s6)^2 + norm(t1-t2)^2 + norm(t3-t4)^2 + norm(t5-t6)^2; +ra = norm(s2)^2 + norm(s4)^2 + norm(s6)^2 + norm(t2)^2 + norm(t4)^2 + norm(t6)^2; + +errs(1) = sqrt(err/ra); +assert(errs(1)