diff --git a/docs/fortran-c.rst b/docs/fortran-c.rst index d6ea3b6..b0049b2 100644 --- a/docs/fortran-c.rst +++ b/docs/fortran-c.rst @@ -330,7 +330,7 @@ denote the Stokeslet given by (x_{3}-y_{3})^2 + \|x-y \|^2 \end{bmatrix} \, , -and $\mathcal{T}^{\textrm{stok}}(x,y)$ denote the Stresslet whose action on +let $\mathcal{T}^{\textrm{stok}}(x,y)$ denote the Stresslet whose action on a vector $v$ is given by .. math:: @@ -343,19 +343,60 @@ a vector $v$ is given by (x_{2}-y_{2})(x_{3}-y_{3}) \\ (x_{3}-y_{3})(x_{1}-y_{1}) & (x_{3}-y_{3})(x_{2}-y_{2}) & (x_{3}-y_{3})^2 - \end{bmatrix} \, . + \end{bmatrix} \, , + +let $\mathcal{R}^{\textrm{stok}}(x,y)$ denote the Rotlet whose action on +a vector $v$ is given by + +.. math:: + v\cdot \mathcal{R}^{\textrm{stok}}(x,y) = + \frac{v \cdot (x-y)}{4\pi\|x-y \|^3} + \begin{bmatrix} + 1 & 0 & 0 \\ + 0 & 1 & 0 \\ + 0 & 0 & 1 + \end{bmatrix} \, , + +and $\mathcal{D}^{\textrm{stok}}(x,y)$ denote the symmetric part of Doublet whose action on +a vector $v$ is given by + +.. math:: + v\cdot \mathcal{D}^{\textrm{stok}}(x,y) = + \frac{3 v \cdot (x-y)}{4\pi\|x-y \|^5} + \begin{bmatrix} + (x_{1}-y_{1})^2 & (x_{1}-y_{1})(x_{2}-y_{2}) & + (x_{1}-y_{1})(x_{3}-y_{3}) \\ + (x_{2}-y_{2})(x_{1}-y_{1}) & (x_{2}-y_{2})^2 & + (x_{2}-y_{2})(x_{3}-y_{3}) \\ + (x_{3}-y_{3})(x_{1}-y_{1}) & (x_{3}-y_{3})(x_{2}-y_{2}) & + (x_{3}-y_{3})^2 + \end{bmatrix} - \\ + \frac{1}{4\pi\|x-y \|^3} + \begin{bmatrix} + v_1(x_{1}-y_{1}) & v_2(x_{1}-y_{1}) & v_3(x_{1}-y_{1}) \\ + v_2(x_{2}-y_{2}) & v_2(x_{2}-y_{2}) & v_3(x_{2}-y_{2}) \\ + v_3(x_{3}-y_{3}) & v_3(x_{3}-y_{3}) & v_3(x_{3}-y_{3}) + \end{bmatrix} \, . The Stokes FMM evaluates the following velocity, its gradient and the associated pressure .. math:: - u(x) = \sum_{m=1}^{N} \mathcal{G}^{\textrm{stok}}(x,x_{j}) \sigma_{j} + \nu_{j} \cdot \mathcal{T}^{\textrm{stok}}(x,x_{j}) \cdot \mu_{j} \, . + u(x) = \sum_{m=1}^{N} \mathcal{G}^{\textrm{stok}}(x,x_{j}) \sigma_{j} + \nu_{j} \cdot \mathcal{T}^{\textrm{stok}}(x,x_{j}) \cdot \mu_{j} + + \nu^{r}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{r}_{j} - \mu^{r}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \nu^{r}_{j} + \\ + \nu^{d}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{d}_{j} - \mu^{d}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \nu^{d}_{j} + + \nu^{d}_{j} \cdot \mathcal{D}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{d}_{j} \, . Here $x_{j}$ are the source locations, $\sigma_{j}$ are the Stokeslet densities, $\nu_{j}$ are the stresslet orientation vectors, $\mu_{j}$ -are the stresslet densities, and the locations $x$ +are the stresslet densities, +$\nu^{r}_{j}$ are the rotlet orientation vectors, $\mu^{r}_{j}$ +are the rotlet densities, +$\nu^{d}_{j}$ are the doublet orientation vectors, $\mu^{d}_{j}$ +are the doublet densities, +and the locations $x$ at which the velocity and its gradient are evaluated are referred to as the evaluation points. @@ -370,7 +411,7 @@ gradients at the evaluation points. .. code:: - subroutine stfmm3d(nd,eps,nsource,source,ifstoklet,stoklet,ifstrslet,strslet,strsvec,ifppreg,pot,pre,grad,ntarg,targ,ifppregtarg,pottarg,pretarg,gradtarg,ier) + subroutine stfmm3d(nd,eps,nsource,source,ifstoklet,stoklet,ifstrslet,strslet,strsvec,ifrotlet,rotlet,rotvec,ifdoublet,doublet,doubvec,ifppreg,pot,pre,grad,ntarg,targ,ifppregtarg,pottarg,pretarg,gradtarg,ier) Input arguments: @@ -396,6 +437,22 @@ Input arguments: Stresslet strengths, $\mu_{j}$ - strsvec: double precision(nd,3,nsource) Stresslet orientation vectors, $\nu_{j}$ + - ifrotlet: integer + Flag for including Rotlet ($\mu^{r}_{j},\nu^{r}_{j}$) term in interaction kernel + Rotlet term will be included if ifrotlet + = 1 + - rotlet: double precision(nd,3,nsource) + Rotlet strengths, $\mu^{r}_{j}$ + - rotvec: double precision(nd,3,nsource) + Rotlet orientation vectors, $\nu^{r}_{j}$ + - ifdoublet: integer + Flag for including Doublet ($\mu^{d}_{j},\nu^{d}_{j}$) term in interaction kernel + Doublet term will be included if ifdoublet + = 1 + - doublet: double precision(nd,3,nsource) + Doublet strengths, $\mu^{d}_{j}$ + - doubvec: double precision(nd,3,nsource) + Doublet orientation vectors, $\nu^{d}_{j}$ - ifppreg: integer | Flag for computing velocity, pressure and/or gradients at source locations | ifppreg = 1, compute velocity diff --git a/docs/genpdfmanual.sh b/docs/genpdfmanual.sh index 9002157..c1235e8 100755 --- a/docs/genpdfmanual.sh +++ b/docs/genpdfmanual.sh @@ -3,4 +3,4 @@ # Barnett 12/6/17 make latexpdf -mv _build/latex/finufft.pdf ../finufft-manual.pdf +mv _build/latex/fmm3d.pdf ../fmm3d_manual.pdf diff --git a/docs/julia.rst b/docs/julia.rst index 521af4f..2655ffe 100644 --- a/docs/julia.rst +++ b/docs/julia.rst @@ -204,7 +204,7 @@ denote the Stokeslet given by (x_{3}-y_{3})^2 + \|x-y \|^2 \end{bmatrix} \, , -and $\mathcal{T}^{\textrm{stok}}(x,y)$ denote the Stresslet whose action on +let $\mathcal{T}^{\textrm{stok}}(x,y)$ denote the Stresslet whose action on a vector $v$ is given by .. math:: @@ -217,7 +217,40 @@ a vector $v$ is given by (x_{2}-y_{2})(x_{3}-y_{3}) \\ (x_{3}-y_{3})(x_{1}-y_{1}) & (x_{3}-y_{3})(x_{2}-y_{2}) & (x_{3}-y_{3})^2 - \end{bmatrix} \, . + \end{bmatrix} \, , + +let $\mathcal{R}^{\textrm{stok}}(x,y)$ denote the Rotlet whose action on +a vector $v$ is given by + +.. math:: + v\cdot \mathcal{R}^{\textrm{stok}}(x,y) = + \frac{v \cdot (x-y)}{4\pi\|x-y \|^3} + \begin{bmatrix} + 1 & 0 & 0 \\ + 0 & 1 & 0 \\ + 0 & 0 & 1 + \end{bmatrix} \, , + +and $\mathcal{D}^{\textrm{stok}}(x,y)$ denote the symmetric part of Doublet whose action on +a vector $v$ is given by + +.. math:: + v\cdot \mathcal{D}^{\textrm{stok}}(x,y) = + \frac{3 v \cdot (x-y)}{4\pi\|x-y \|^5} + \begin{bmatrix} + (x_{1}-y_{1})^2 & (x_{1}-y_{1})(x_{2}-y_{2}) & + (x_{1}-y_{1})(x_{3}-y_{3}) \\ + (x_{2}-y_{2})(x_{1}-y_{1}) & (x_{2}-y_{2})^2 & + (x_{2}-y_{2})(x_{3}-y_{3}) \\ + (x_{3}-y_{3})(x_{1}-y_{1}) & (x_{3}-y_{3})(x_{2}-y_{2}) & + (x_{3}-y_{3})^2 + \end{bmatrix} - \\ + \frac{1}{4\pi\|x-y \|^3} + \begin{bmatrix} + v_1(x_{1}-y_{1}) & v_2(x_{1}-y_{1}) & v_3(x_{1}-y_{1}) \\ + v_2(x_{2}-y_{2}) & v_2(x_{2}-y_{2}) & v_3(x_{2}-y_{2}) \\ + v_3(x_{3}-y_{3}) & v_3(x_{3}-y_{3}) & v_3(x_{3}-y_{3}) + \end{bmatrix} \, . This subroutine computes the N-body Stokes interactions, its gradients and the corresponding pressure @@ -225,11 +258,19 @@ in three dimensions given by .. math:: - u(x) = \sum_{m=1}^{N} \mathcal{G}^{\textrm{stok}}(x,x_{j}) \sigma_{j} + \nu_{j} \cdot \mathcal{T}^{\textrm{stok}}(x,x_{j}) \cdot \mu_{j} + u(x) = \sum_{m=1}^{N} \mathcal{G}^{\textrm{stok}}(x,x_{j}) \sigma_{j} + \nu_{j} \cdot \mathcal{T}^{\textrm{stok}}(x,x_{j}) \cdot \mu_{j} + + \nu^{r}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{r}_{j} - \mu^{r}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \nu^{r}_{j} + \\ + \nu^{d}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{d}_{j} - \mu^{d}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \nu^{d}_{j} + + \nu^{d}_{j} \cdot \mathcal{D}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{d}_{j} \, . + where $\sigma_{j}$ are the Stokeslet densities, $\nu_{j}$ are the stresslet orientation vectors, $\mu_{j}$ -are the stresslet densities, and +are the stresslet densities, +$\nu^{r}_{j}$ are the rotlet orientation vectors, $\mu^{r}_{j}$ +are the rotlet densities, +$\nu^{d}_{j}$ are the doublet orientation vectors, $\mu^{d}_{j}$ +are the doublet densities, and $x_{j}$ are the source locations. When $x=x_{j}$, the term corresponding to $x_{j}$ is dropped from the sum. @@ -237,7 +278,8 @@ from the sum. .. code:: julia vals = stfmm3d(eps,sources;stoklet=nothing,strslet=nothing, - strsvec=nothing,targets=nothing,ppreg=0, + strsvec=nothing,rotlet=nothing,rotvec=nothing, + doublet=nothing,doubvec=nothing,targets=nothing,ppreg=0, ppregt=0,nd=1) Wrapper for fast multipole implementation for Stokes N-body @@ -255,6 +297,14 @@ Args: stresslet strengths ($mu_{j}$ above) - strsvec: float(nd,3,n) or float(3,n) stresslet orientations ($nu_{j}$ above) +- rotlet: float(nd,3,n) or float(3,n) + rotlet strengths ($mu^{r}_{j}$ above) +- rotvec: float(nd,3,n) or float(3,n) + rotlet orientations ($nu^{r}_{j}$ above) +- doublet: float(nd,3,n) or float(3,n) + doublet strengths ($mu^{d}_{j}$ above) +- doubvec: float(nd,3,n) or float(3,n) + doublet orientations ($nu^{d}_{j}$ above) - targets: float(3,nt) target locations (x) - ifppreg: integer @@ -286,7 +336,8 @@ target locations. .. code:: julia vals = st3ddir(sources,targets;stoklet=nothing,strslet=nothing, - strsvec=nothing,ppregt=0,nd=1,thresh=1e-16) + strsvec=nothing,rotlet=nothing,rotvec=nothing, + doublet=nothing,doubvec=nothing,ppregt=0,nd=1,thresh=1e-16) ------------------------------------------------------------------ diff --git a/docs/matlab.rst b/docs/matlab.rst index 45d9782..3dc31f7 100644 --- a/docs/matlab.rst +++ b/docs/matlab.rst @@ -205,7 +205,7 @@ denote the Stokeslet given by (x_{3}-y_{3})^2 + \|x-y \|^2 \end{bmatrix} \, , -and $\mathcal{T}^{\textrm{stok}}(x,y)$ denote the Stresslet whose action on +let $\mathcal{T}^{\textrm{stok}}(x,y)$ denote the Stresslet whose action on a vector $v$ is given by .. math:: @@ -218,19 +218,59 @@ a vector $v$ is given by (x_{2}-y_{2})(x_{3}-y_{3}) \\ (x_{3}-y_{3})(x_{1}-y_{1}) & (x_{3}-y_{3})(x_{2}-y_{2}) & (x_{3}-y_{3})^2 - \end{bmatrix} \, . + \end{bmatrix} \, , + +let $\mathcal{R}^{\textrm{stok}}(x,y)$ denote the Rotlet whose action on +a vector $v$ is given by + +.. math:: + v\cdot \mathcal{R}^{\textrm{stok}}(x,y) = + \frac{v \cdot (x-y)}{4\pi\|x-y \|^3} + \begin{bmatrix} + 1 & 0 & 0 \\ + 0 & 1 & 0 \\ + 0 & 0 & 1 + \end{bmatrix} \, , + +and $\mathcal{D}^{\textrm{stok}}(x,y)$ denote the symmetric part of Doublet whose action on +a vector $v$ is given by + +.. math:: + v\cdot \mathcal{D}^{\textrm{stok}}(x,y) = + \frac{3 v \cdot (x-y)}{4\pi\|x-y \|^5} + \begin{bmatrix} + (x_{1}-y_{1})^2 & (x_{1}-y_{1})(x_{2}-y_{2}) & + (x_{1}-y_{1})(x_{3}-y_{3}) \\ + (x_{2}-y_{2})(x_{1}-y_{1}) & (x_{2}-y_{2})^2 & + (x_{2}-y_{2})(x_{3}-y_{3}) \\ + (x_{3}-y_{3})(x_{1}-y_{1}) & (x_{3}-y_{3})(x_{2}-y_{2}) & + (x_{3}-y_{3})^2 + \end{bmatrix} - \\ + \frac{1}{4\pi\|x-y \|^3} + \begin{bmatrix} + v_1(x_{1}-y_{1}) & v_2(x_{1}-y_{1}) & v_3(x_{1}-y_{1}) \\ + v_2(x_{2}-y_{2}) & v_2(x_{2}-y_{2}) & v_3(x_{2}-y_{2}) \\ + v_3(x_{3}-y_{3}) & v_3(x_{3}-y_{3}) & v_3(x_{3}-y_{3}) + \end{bmatrix} \, . This subroutine computes the N-body Stokes -interactions, its gradients and the corresponding pressure -in three dimensions given by - +interactions, its gradients and the corresponding pressure +in three dimensions given by + .. math:: - u(x) = \sum_{m=1}^{N} \mathcal{G}^{\textrm{stok}}(x,x_{j}) \sigma_{j} + \nu_{j} \cdot \mathcal{T}^{\textrm{stok}}(x,x_{j}) \cdot \mu_{j} + u(x) = \sum_{m=1}^{N} \mathcal{G}^{\textrm{stok}}(x,x_{j}) \sigma_{j} + \nu_{j} \cdot \mathcal{T}^{\textrm{stok}}(x,x_{j}) \cdot \mu_{j} + + \nu^{r}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{r}_{j} - \mu^{r}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \nu^{r}_{j} + \\ + \nu^{d}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{d}_{j} - \mu^{d}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \nu^{d}_{j} + + \nu^{d}_{j} \cdot \mathcal{D}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{d}_{j} \, . where $\sigma_{j}$ are the Stokeslet densities, $\nu_{j}$ are the stresslet orientation vectors, $\mu_{j}$ -are the stresslet densities, and +are the stresslet densities, +$\nu^{r}_{j}$ are the rotlet orientation vectors, $\mu^{r}_{j}$ +are the rotlet densities, +$\nu^{d}_{j}$ are the doublet orientation vectors, $\mu^{d}_{j}$ +are the doublet densities, and $x_{j}$ are the source locations. When $x=x_{j}$, the term corresponding to $x_{j}$ is dropped from the sum. @@ -263,6 +303,19 @@ Args: * srcinfo.strsvec: double(nd,3,n) Stresslet orientiation vectors, $\nu_{j}$ (optional default - term corresponding to stresslet dropped) + * srcinfo.rotlet: double(nd,3,n) + Rotlet densities, $\mu^{r}_{j}$ (optional + default - term corresponding to rotlet dropped) + * srcinfo.rotvec: double(nd,3,n) + Rotlet orientiation vectors, $\nu^{r}_{j}$ (optional + default - term corresponding to rotlet dropped) + * srcinfo.doublet: double(nd,3,n) + Doublet densities, $\mu^{d}_{j}$ (optional + default - term corresponding to doublet dropped) + * srcinfo.doubvec: double(nd,3,n) + Doublet orientiation vectors, $\nu^{d}_{j}$ (optional + default - term corresponding to doublet dropped) + - ifppreg: integer | source eval flag diff --git a/docs/python.rst b/docs/python.rst index 6feee26..717cb1c 100644 --- a/docs/python.rst +++ b/docs/python.rst @@ -194,7 +194,7 @@ denote the Stokeslet given by (x_{3}-y_{3})^2 + \|x-y \|^2 \end{bmatrix} \, , -and $\mathcal{T}^{\textrm{stok}}(x,y)$ denote the Stresslet whose action on +let $\mathcal{T}^{\textrm{stok}}(x,y)$ denote the Stresslet whose action on a vector $v$ is given by .. math:: @@ -207,26 +207,67 @@ a vector $v$ is given by (x_{2}-y_{2})(x_{3}-y_{3}) \\ (x_{3}-y_{3})(x_{1}-y_{1}) & (x_{3}-y_{3})(x_{2}-y_{2}) & (x_{3}-y_{3})^2 - \end{bmatrix} \, . + \end{bmatrix} \, , + +let $\mathcal{R}^{\textrm{stok}}(x,y)$ denote the Rotlet whose action on +a vector $v$ is given by + +.. math:: + v\cdot \mathcal{R}^{\textrm{stok}}(x,y) = + \frac{v \cdot (x-y)}{4\pi\|x-y \|^3} + \begin{bmatrix} + 1 & 0 & 0 \\ + 0 & 1 & 0 \\ + 0 & 0 & 1 + \end{bmatrix} \, , + +and $\mathcal{D}^{\textrm{stok}}(x,y)$ denote the symmetric part of Doublet whose action on +a vector $v$ is given by + +.. math:: + v\cdot \mathcal{D}^{\textrm{stok}}(x,y) = + \frac{3 v \cdot (x-y)}{4\pi\|x-y \|^5} + \begin{bmatrix} + (x_{1}-y_{1})^2 & (x_{1}-y_{1})(x_{2}-y_{2}) & + (x_{1}-y_{1})(x_{3}-y_{3}) \\ + (x_{2}-y_{2})(x_{1}-y_{1}) & (x_{2}-y_{2})^2 & + (x_{2}-y_{2})(x_{3}-y_{3}) \\ + (x_{3}-y_{3})(x_{1}-y_{1}) & (x_{3}-y_{3})(x_{2}-y_{2}) & + (x_{3}-y_{3})^2 + \end{bmatrix} - \\ + \frac{1}{4\pi\|x-y \|^3} + \begin{bmatrix} + v_1(x_{1}-y_{1}) & v_2(x_{1}-y_{1}) & v_3(x_{1}-y_{1}) \\ + v_2(x_{2}-y_{2}) & v_2(x_{2}-y_{2}) & v_3(x_{2}-y_{2}) \\ + v_3(x_{3}-y_{3}) & v_3(x_{3}-y_{3}) & v_3(x_{3}-y_{3}) + \end{bmatrix} \, . This subroutine computes the N-body Stokes interactions, its gradients and the corresponding pressure in three dimensions given by - + + .. math:: - u(x) = \sum_{m=1}^{N} \mathcal{G}^{\textrm{stok}}(x,x_{j}) \sigma_{j} + \nu_{j} \cdot \mathcal{T}^{\textrm{stok}}(x,x_{j}) \cdot \mu_{j} + u(x) = \sum_{m=1}^{N} \mathcal{G}^{\textrm{stok}}(x,x_{j}) \sigma_{j} + \nu_{j} \cdot \mathcal{T}^{\textrm{stok}}(x,x_{j}) \cdot \mu_{j} + + \nu^{r}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{r}_{j} - \mu^{r}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \nu^{r}_{j} + \\ + \nu^{d}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{d}_{j} - \mu^{d}_{j} \cdot \mathcal{R}^{\textrm{stok}}(x,x_{j}) \cdot \nu^{d}_{j} + + \nu^{d}_{j} \cdot \mathcal{D}^{\textrm{stok}}(x,x_{j}) \cdot \mu^{d}_{j} \, . where $\sigma_{j}$ are the Stokeslet densities, $\nu_{j}$ are the stresslet orientation vectors, $\mu_{j}$ -are the stresslet densities, and +are the stresslet densities, +$\nu^{r}_{j}$ are the rotlet orientation vectors, $\mu^{r}_{j}$ +are the rotlet densities, +$\nu^{d}_{j}$ are the doublet orientation vectors, $\mu^{d}_{j}$ +are the doublet densities, and $x_{j}$ are the source locations. When $x=x_{j}$, the term corresponding to $x_{j}$ is dropped from the sum. .. code:: python - def stfmm3d(*,eps,sources,stoklet=None,strslet=None,strsvec=None,targets=None,ifppreg=0,ifppregtarg=0,nd=1) + def stfmm3d(*,eps,sources,stoklet=None,strslet=None,strsvec=None,rotlet=None,rotvec=None,doublet=None,doubvec=None,targets=None,ifppreg=0,ifppregtarg=0,nd=1) Wrapper for fast multipole implementation for Stokes N-body interactions. @@ -243,6 +284,14 @@ Args: stresslet strengths ($\mu_{j}$ above) - strsvec: float(nd,3,n) or float(3,n) stresslet orientations ($\nu_{j}$ above) +- rotlet: float(nd,3,n) or float(3,n) + rotlet strengths ($\mu^{r}_{j}$ above) +- rotvec: float(nd,3,n) or float(3,n) + rotlet orientations ($\nu^{r}_{j}$ above) +- doublet: float(nd,3,n) or float(3,n) + doublet strengths ($\mu^{d}_{j}$ above) +- doubvec: float(nd,3,n) or float(3,n) + doublet orientations ($\nu^{d}_{j}$ above) - targets: float(3,nt) target locations (x) - ifppreg: integer @@ -273,7 +322,7 @@ target locations. .. code:: python - def st3ddir(*,eps,sources,stoklet=None,strslet=None,strsvec=None,targets=None,ifppreg=0,ifppregtarg=0,nd=1,thresh=1e-16): + def st3ddir(*,eps,sources,stoklet=None,strslet=None,strsvec=None,rotlet=None,rotvec=None,doublet=None,doubvec=None,targets=None,ifppreg=0,ifppregtarg=0,nd=1,thresh=1e-16): ------------------------------------------------------------------ diff --git a/examples/stfmm3d_example.f b/examples/stfmm3d_example.f index 570caa3..c8400b1 100644 --- a/examples/stfmm3d_example.f +++ b/examples/stfmm3d_example.f @@ -13,6 +13,8 @@ real *8, allocatable :: stoklet(:,:), strslet(:,:) real *8, allocatable :: strsvec(:,:) + real *8, allocatable :: rotlet(:,:), rotvec(:,:) + real *8, allocatable :: doublet(:,:), doubvec(:,:) integer *8 ipass(10) complex *16 eye @@ -30,6 +32,8 @@ allocate(source(3,ns),targ(3,nt)) allocate(stoklet(3,ns),strslet(3,ns),strsvec(3,ns)) + allocate(rotlet(3,ns),rotvec(3,ns)) + allocate(doublet(3,ns),doubvec(3,ns)) c @@ -67,14 +71,17 @@ eps = 1d-9 ifstoklet = 1 ifstrslet = 1 + ifrotlet = 0 + ifdoublet = 0 ifppreg = 2 ifppregtarg = 0 call cpu_time(t1) c$ t1 = omp_get_wtime() call stfmm3d(nd,eps,ns,source,ifstoklet,stoklet, - 1 ifstrslet,strslet,strsvec,ifppreg,pot,pre,grad, - 2 nt,targ,ifppregtarg,pottarg,pretarg,gradtarg,ier) + 1 ifstrslet,strslet,strsvec,ifrotlet,rotlet,rotvec, + 2 ifdoublet,doublet,doubvec,ifppreg,pot,pre,grad, + 3 nt,targ,ifppregtarg,pottarg,pretarg,gradtarg,ier) call cpu_time(t2) c$ t2 = omp_get_wtime() diff --git a/julia/src/FMM3D.jl b/julia/src/FMM3D.jl index 29afc28..24d12be 100644 --- a/julia/src/FMM3D.jl +++ b/julia/src/FMM3D.jl @@ -1485,7 +1485,7 @@ function stfmm3d(eps::Float64, # $ ifppregtarg, pottarg, pretarg, gradtarg,ier) - ccall((:stfmm3d_new_,libfmm3d),Cvoid,(Fi,Fd,Fi,Fd,Fi,Fd, + ccall((:stfmm3d_,libfmm3d),Cvoid,(Fi,Fd,Fi,Fd,Fi,Fd, Fi,Fd,Fd,Fi,Fd,Fd,Fi,Fd,Fd, Fi,Fd,Fd,Fd,Fi, Fd,Fi,Fd,Fd,Fd,Fi), diff --git a/matlab/fmm3d.c b/matlab/fmm3d.c index f46d97a..d9bebfa 100644 --- a/matlab/fmm3d.c +++ b/matlab/fmm3d.c @@ -1090,7 +1090,7 @@ mxWrapReturnZDef_single (mxWrapReturn_single_dcomplex, dcomplex, #define MWF77_l3ddirectcdh L3DDIRECTCDH #define MWF77_emfmm3d EMFMM3D #define MWF77_em3ddirect EM3DDIRECT -#define MWF77_stfmm3d_new STFMM3D_NEW +#define MWF77_stfmm3d STFMM3D #define MWF77_st3ddirectstokstrsrotdoubg ST3DDIRECTSTOKSTRSROTDOUBG #elif defined(MWF77_UNDERSCORE1) #define MWF77_hndiv hndiv_ @@ -1114,7 +1114,7 @@ mxWrapReturnZDef_single (mxWrapReturn_single_dcomplex, dcomplex, #define MWF77_l3ddirectcdh l3ddirectcdh_ #define MWF77_emfmm3d emfmm3d_ #define MWF77_em3ddirect em3ddirect_ -#define MWF77_stfmm3d_new stfmm3d_new_ +#define MWF77_stfmm3d stfmm3d_ #define MWF77_st3ddirectstokstrsrotdoubg st3ddirectstokstrsrotdoubg_ #elif defined(MWF77_UNDERSCORE0) #define MWF77_hndiv hndiv @@ -1138,7 +1138,7 @@ mxWrapReturnZDef_single (mxWrapReturn_single_dcomplex, dcomplex, #define MWF77_l3ddirectcdh l3ddirectcdh #define MWF77_emfmm3d emfmm3d #define MWF77_em3ddirect em3ddirect -#define MWF77_stfmm3d_new stfmm3d_new +#define MWF77_stfmm3d stfmm3d #define MWF77_st3ddirectstokstrsrotdoubg st3ddirectstokstrsrotdoubg #else /* f2c convention */ #define MWF77_hndiv hndiv_ @@ -1162,7 +1162,7 @@ mxWrapReturnZDef_single (mxWrapReturn_single_dcomplex, dcomplex, #define MWF77_l3ddirectcdh l3ddirectcdh_ #define MWF77_emfmm3d emfmm3d_ #define MWF77_em3ddirect em3ddirect_ -#define MWF77_stfmm3d_new stfmm3d_new__ +#define MWF77_stfmm3d stfmm3d_ #define MWF77_st3ddirectstokstrsrotdoubg st3ddirectstokstrsrotdoubg_ #endif @@ -1195,7 +1195,7 @@ MWF77_RETURN MWF77_l3ddirectdh(int64_t*, double*, double*, int64_t*, double*, in MWF77_RETURN MWF77_l3ddirectcdh(int64_t*, double*, double*, double*, int64_t*, double*, int64_t*, double*, double*, double*, double*); MWF77_RETURN MWF77_emfmm3d(int64_t*, double*, dcomplex*, int64_t*, double*, int64_t*, dcomplex*, int64_t*, dcomplex*, int64_t*, dcomplex*, int64_t*, double*, int64_t*, dcomplex*, int64_t*, dcomplex*, int64_t*, dcomplex*, int64_t*); MWF77_RETURN MWF77_em3ddirect(int64_t*, dcomplex*, int64_t*, double*, int64_t*, dcomplex*, int64_t*, dcomplex*, int64_t*, dcomplex*, int64_t*, double*, int64_t*, dcomplex*, int64_t*, dcomplex*, int64_t*, dcomplex*, double*); -MWF77_RETURN MWF77_stfmm3d_new(int64_t*, double*, int64_t*, double*, int64_t*, double*, int64_t*, double*, double*, int64_t*, double*, double*, int64_t*, double*, double*, int64_t*, double*, double*, double*, int64_t*, double*, int64_t*, double*, double*, double*, int64_t*); +MWF77_RETURN MWF77_stfmm3d(int64_t*, double*, int64_t*, double*, int64_t*, double*, int64_t*, double*, double*, int64_t*, double*, double*, int64_t*, double*, double*, int64_t*, double*, double*, double*, int64_t*, double*, int64_t*, double*, double*, double*, int64_t*); MWF77_RETURN MWF77_st3ddirectstokstrsrotdoubg(int64_t*, double*, double*, int64_t*, double*, double*, int64_t*, double*, double*, int64_t*, double*, double*, int64_t*, double*, int64_t*, double*, double*, double*, double*); #ifdef __cplusplus @@ -6099,9 +6099,9 @@ void mexStub21(int nlhs, mxArray* plhs[], } /* ---- fmm3d.mw: 1134 ---- - * stfmm3d_new(int64_t[1] nd, double[1] eps, int64_t[1] ns, double[3, ns] sources, int64_t[1] ifstoklet, double[nd3, ns_stok] stoklet, int64_t[1] ifstrslet, double[nd3, ns_strs] strslet, double[nd3, ns_strs] strsvec, int64_t[1] ifrotlet, double[nd3, ns_rot] rotlet, double[nd3, ns_rot] rotvec, int64_t[1] ifdoublet, double[nd3, ns_doub] doublet, double[nd3, ns_doub] doubvec, int64_t[1] ifppreg, inout double[nd3, ns_pot] pot, inout double[nd, ns_pre] pre, inout double[nd9, ns_grad] grad, int64_t[1] nt, double[3, nt] targ, int64_t[1] ifppregtarg, inout double[nd3, nt_pot] pottarg, inout double[nd, nt_pre] pretarg, inout double[nd9, nt_grad] gradtarg, inout int64_t[1] ier); + * stfmm3d(int64_t[1] nd, double[1] eps, int64_t[1] ns, double[3, ns] sources, int64_t[1] ifstoklet, double[nd3, ns_stok] stoklet, int64_t[1] ifstrslet, double[nd3, ns_strs] strslet, double[nd3, ns_strs] strsvec, int64_t[1] ifrotlet, double[nd3, ns_rot] rotlet, double[nd3, ns_rot] rotvec, int64_t[1] ifdoublet, double[nd3, ns_doub] doublet, double[nd3, ns_doub] doubvec, int64_t[1] ifppreg, inout double[nd3, ns_pot] pot, inout double[nd, ns_pre] pre, inout double[nd9, ns_grad] grad, int64_t[1] nt, double[3, nt] targ, int64_t[1] ifppregtarg, inout double[nd3, nt_pot] pottarg, inout double[nd, nt_pre] pretarg, inout double[nd9, nt_grad] gradtarg, inout int64_t[1] ier); */ -static const char* stubids22_ = "stfmm3d_new(i int64_t[x], i double[x], i int64_t[x], i double[xx], i int64_t[x], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], io double[xx], io double[xx], io double[xx], i int64_t[x], i double[xx], i int64_t[x], io double[xx], io double[xx], io double[xx], io int64_t[x])"; +static const char* stubids22_ = "stfmm3d(i int64_t[x], i double[x], i int64_t[x], i double[xx], i int64_t[x], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], io double[xx], io double[xx], io double[xx], i int64_t[x], i double[xx], i int64_t[x], io double[xx], io double[xx], io double[xx], io int64_t[x])"; void mexStub22(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) @@ -6559,7 +6559,7 @@ void mexStub22(int nlhs, mxArray* plhs[], in25_ = NULL; if (mexprofrecord_) mexprofrecord_[22]++; - MWF77_stfmm3d_new(in0_, in1_, in2_, in3_, in4_, in5_, in6_, in7_, in8_, in9_, in10_, in11_, in12_, in13_, in14_, in15_, in16_, in17_, in18_, in19_, in20_, in21_, in22_, in23_, in24_, in25_); + MWF77_stfmm3d(in0_, in1_, in2_, in3_, in4_, in5_, in6_, in7_, in8_, in9_, in10_, in11_, in12_, in13_, in14_, in15_, in16_, in17_, in18_, in19_, in20_, in21_, in22_, in23_, in24_, in25_); plhs[0] = mxCreateDoubleMatrix(dim50_, dim51_, mxREAL); mxWrapCopy_double(plhs[0], in16_, dim50_*dim51_); plhs[1] = mxCreateDoubleMatrix(dim52_, dim53_, mxREAL); diff --git a/matlab/fmm3d.mw b/matlab/fmm3d.mw index 6135f15..82c0394 100644 --- a/matlab/fmm3d.mw +++ b/matlab/fmm3d.mw @@ -1131,7 +1131,7 @@ end nd9 = 9*nd; ier = 0; - # FORTRAN stfmm3d_new(int[1] nd, double[1] eps, int[1] ns, double[3,ns] sources, int[1] ifstoklet, double[nd3,ns_stok] stoklet, int[1] ifstrslet, double[nd3,ns_strs] strslet, double[nd3,ns_strs] strsvec, int[1] ifrotlet, double[nd3,ns_rot] rotlet, double[nd3,ns_rot] rotvec, int[1] ifdoublet, double[nd3,ns_doub] doublet, double[nd3,ns_doub] doubvec, int[1] ifppreg, inout double[nd3,ns_pot] pot, inout double[nd,ns_pre] pre, inout double[nd9,ns_grad] grad, int[1] nt, double[3,nt] targ, int[1] ifppregtarg, inout double[nd3,nt_pot] pottarg, inout double [nd,nt_pre] pretarg, inout double[nd9,nt_grad] gradtarg, inout int[1] ier); + # FORTRAN stfmm3d(int[1] nd, double[1] eps, int[1] ns, double[3,ns] sources, int[1] ifstoklet, double[nd3,ns_stok] stoklet, int[1] ifstrslet, double[nd3,ns_strs] strslet, double[nd3,ns_strs] strsvec, int[1] ifrotlet, double[nd3,ns_rot] rotlet, double[nd3,ns_rot] rotvec, int[1] ifdoublet, double[nd3,ns_doub] doublet, double[nd3,ns_doub] doubvec, int[1] ifppreg, inout double[nd3,ns_pot] pot, inout double[nd,ns_pre] pre, inout double[nd9,ns_grad] grad, int[1] nt, double[3,nt] targ, int[1] ifppregtarg, inout double[nd3,nt_pot] pottarg, inout double [nd,nt_pre] pretarg, inout double[nd9,nt_grad] gradtarg, inout int[1] ier); U.pot = []; U.pre = []; diff --git a/matlab/stfmm3d.m b/matlab/stfmm3d.m index 486fd2f..32e52ca 100644 --- a/matlab/stfmm3d.m +++ b/matlab/stfmm3d.m @@ -231,7 +231,7 @@ nd9 = 9*nd; ier = 0; - mex_id_ = 'stfmm3d_new(i int64_t[x], i double[x], i int64_t[x], i double[xx], i int64_t[x], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], io double[xx], io double[xx], io double[xx], i int64_t[x], i double[xx], i int64_t[x], io double[xx], io double[xx], io double[xx], io int64_t[x])'; + mex_id_ = 'stfmm3d(i int64_t[x], i double[x], i int64_t[x], i double[xx], i int64_t[x], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], i double[xx], i double[xx], i int64_t[x], io double[xx], io double[xx], io double[xx], i int64_t[x], i double[xx], i int64_t[x], io double[xx], io double[xx], io double[xx], io int64_t[x])'; [pot, pre, grad, pottarg, pretarg, gradtarg, ier] = fmm3d(mex_id_, nd, eps, ns, sources, ifstoklet, stoklet, ifstrslet, strslet, strsvec, ifrotlet, rotlet, rotvec, ifdoublet, doublet, doubvec, ifppreg, pot, pre, grad, nt, targ, ifppregtarg, pottarg, pretarg, gradtarg, ier, 1, 1, 1, 3, ns, 1, nd3, ns_stok, 1, nd3, ns_strs, nd3, ns_strs, 1, nd3, ns_rot, nd3, ns_rot, 1, nd3, ns_doub, nd3, ns_doub, 1, nd3, ns_pot, nd, ns_pre, nd9, ns_grad, 1, 3, nt, 1, nd3, nt_pot, nd, nt_pre, nd9, nt_grad, 1); U.pot = []; diff --git a/python/fmm3dpy/fmm3d.py b/python/fmm3dpy/fmm3d.py index d096d75..7890fa6 100644 --- a/python/fmm3dpy/fmm3d.py +++ b/python/fmm3dpy/fmm3d.py @@ -748,7 +748,7 @@ def stfmm3d(*,eps,sources,stoklet=None,strslet=None,strsvec=None,rotlet=None,rot doublet = np.zeros([nd,3,ns],dtype='double') doubvec = np.zeros([nd,3,ns],dtype='double') - out.pot,out.pre,out.grad,out.pottarg,out.pretarg,out.gradtarg,out.ier = stfmm.stfmm3d_new(eps,sources,ifstoklet,stoklet,ifstrslet,strslet,strsvec,ifrotlet,rotlet,rotvec,ifdoublet,doublet,doubvec,ifppreg,targets,ifppregtarg,nd,ns,nt) + out.pot,out.pre,out.grad,out.pottarg,out.pretarg,out.gradtarg,out.ier = stfmm.stfmm3d(eps,sources,ifstoklet,stoklet,ifstrslet,strslet,strsvec,ifrotlet,rotlet,rotvec,ifdoublet,doublet,doubvec,ifppreg,targets,ifppregtarg,nd,ns,nt) if(ifppreg < 3): out.grad = None diff --git a/src/Stokes/stfmm3d.f b/src/Stokes/stfmm3d.f index 4d2bf7b..604ba52 100644 --- a/src/Stokes/stfmm3d.f +++ b/src/Stokes/stfmm3d.f @@ -2,7 +2,7 @@ c This file contains the Stokes FMM wrappers c c********************************************************************** -c +c c We take the following conventions for the Stokes kernels. c c 1) The dynamic viscosity (mu) is assumed to be 1. @@ -19,475 +19,23 @@ c c The (Type I) stresslet, T_{ijk}, and its associated pressure c tensor, PI_{jk}, we define as -c +c c T_{ijk}(x,y) = 3/(4\pi) r_i r_j r_k/ r^5 -c PI_{jk}(x,y) = -1/(2\pi) delta_{jk}/r^3 + 3/(2\pi) r_j r_k/r^5 - +c PI_{jk}(x,y) = -1/(2\pi) delta_{jk}/r^3 + 3/(2\pi) r_j r_k/r^5 - subroutine stfmm3d(nd, eps, - $ nsource, source, - $ ifstoklet, stoklet, ifstrslet, strslet, strsvec, - $ ifppreg, pot, pre, grad, ntarg, targ, - $ ifppregtarg, pottarg, pretarg, gradtarg,ier) -cf2py intent(in) nd,eps -cf2py intent(in) nsource,source -cf2py intent(in) ifstoklet,stoklet -cf2py intent(in) ifstrslet,strslet,strsvec -cf2py intent(in) ifppreg,ifppregtarg -cf2py intent(in) ntarg,targ -cf2py intent(out) pot,pre,grad -cf2py intent(out) pottarg,pretarg,gradtarg -cf2py intent(out) ier -c -c Stokes FMM in R^{3}: evaluate all pairwise particle -c interactions (ignoring self-interactions) and -c interactions with targs. -c -c This routine computes the sum for the velocity vector, -c -c u_i(x) = sum_m G_{ij}(x,y^{(m)}) sigma^{(m)}_j -c + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k -c -c for x at all of the target locations, and i=1,2,3, -c where sigma^{(m)} is the Stokeslet force, mu^{(m)} is the -c stresslet strength, and nu^{(m)} is the stresslet orientation -c (note that each of these is a 3 vector per source point y^{(m)}). -c Repeated indices are taken as summed over 1,2,3, ie, Einstein -c convention. For x a source point, the self-interaction in the -c sum is omitted. -c -c Optionally, the associated pressure p(x) and 3x3 gradient tensor -c grad u(x) are returned -c -c p(x) = sum_m P_j(x,y^m) sigma^{(m)}_j -c + sum_m PI_{jk}(x,y{(m)}) mu^{(m)}_j nu^{(m)}_k -c -c grad_l u_i(x) = grad_l [sum_m G_{ij}(x,y^m) sigma^{(m)}_j -c + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k] -c -c Note that these two may be combined to get the stress tensor. -c -c----------------------------------------------------------------------- -c INPUT PARAMETERS: -c -c nd: in: integer *8 -c number of densities -c -c eps: in: double precision -c requested precision -c -c nsource in: integer *8 -c number of sources -c -c source in: double precision (3,nsource) -c source(k,j) is the kth component of the jth -c source locations -c -c ifstoklet in: integer *8 -c Stokeslet charge computation flag -c ifstoklet = 1 => include Stokeslet contribution -c otherwise do not -c -c stoklet in: double precision (nd,3,nsource) -c Stokeslet charge strengths (sigma vectors above) -c -c ifstrslet in: integer *8 -c stresslet computation flag -c ifstrslet = 1 => include standard stresslet -c (type I) +c The rotlet, R_{ijk}, and its associated pressure tensor, Q_{jk}, are c -c NOT YET IMPLEMENTED -c -c ifstrslet = 2 => include symmetric stresslet -c (type II) -c ifstrslet = 3 => include rotlet -c ifstrslet = 4 => include Stokes doublet -c otherwise do not include +c R_{ijk}(x,y) = -\delta_{ik} r_j/(4\pi r^3) + \delta_{ij} r_k/(4\pi r^3) +c Q_{jk} = 0; c -c strslet in: double precision (nd,3,nsource) -c stresslet strengths (mu vectors above) -c -c strsvec in: double precision (nd,3,nsource) -c stresslet orientations (nu vectors above) +c The doublet, D_{ijk}, and its associated pressure tensor, L_{jk}, are c -c ifppreg in: integer *8 -c flag for evaluating potential, gradient, and pressure -c at the sources -c ifppreg = 1, only potential -c ifppreg = 2, potential and pressure -c GRADIENT NOT IMPLEMENTED -c ifppreg = 3, potential, pressure, and gradient -c -c ntarg in: integer *8 -c number of targs -c -c targ in: double precision (3,ntarg) -c targ(k,j) is the kth component of the jth -c targ location -c -c ifppregtarg in: integer *8 -c flag for evaluating potential, gradient, and pressure -c at the targets -c ifppregtarg = 1, only potential -c ifppregtarg = 2, potential and pressure -c ifppregtarg = 3, potential, pressure, and gradient -c -c----------------------------------------------------------------------- -c -c OUTPUT parameters: -c -c pot out: double precision(nd,3,nsource) -c velocity at the source locations -c -c pre out: double precision(nd,nsource) -c pressure at the source locations -c -c GRADIENT NOT IMPLEMENTED -c grad out: double precision(nd,3,3,nsource) -c gradient of velocity at the source locations -c grad(l,i,j,k) is the ith component of the -c gradient of the jth component of the velocity -c for the lth density at the kth source location -c -c pottarg out: double precision(nd,3,ntarg) -c velocity at the targets -c -c pretarg out: double precision(nd,ntarg) -c pressure at the targets -c -c gradtarg out: double precision(nd,3,3,ntarg) -c gradient of velocity at the targets -c gradtarg(l,i,j,k) is the ith component of the -c gradient of the jth component of the velocity -c for the lth density at the kth target -c ier out: integer *8 -c error flag -c -c TODO: implement other stresslet options and gradient -c------------------------------------------------------------------ - implicit none - integer *8 nd, ifstoklet, ifstrslet, ntarg - double precision eps - integer *8 nsource, ifppreg, ifppregtarg - double precision source(3, nsource), targ(3, ntarg) - double precision stoklet(nd, 3, nsource), strslet(nd, 3, nsource) - double precision strsvec(nd, 3, nsource) - double precision pot(nd, 3, nsource), pre(nd,nsource) - double precision grad(nd, 3, 3, nsource) - double precision pottarg(nd, 3, ntarg), pretarg(nd,ntarg), - 1 gradtarg(nd, 3, 3, ntarg) - -c local - double precision, allocatable :: charge(:,:,:), dipvec(:,:,:,:) - double precision, allocatable :: potl(:,:,:), gradl(:,:,:,:), - 1 hessl(:,:,:,:), pottargl(:,:,:), gradtargl(:,:,:,:), - 2 hesstargl(:,:,:,:) - double precision :: pt(3), gl(3), hl(6), vel(3), velgrad(3,3) - double precision :: press, pl, pv, dmu(3), dnu(3), sigma(3) - - integer *8 ndl, ifchargel, ifdipolel, ifpghl, ifpghtargl - integer *8 ndper - - integer *8 i, j, ii, ifppreg1, l, npt, ier,iper - - - ndper = 0 - - if (ifstrslet .eq. 1 .or. ifstoklet .eq. 1) then - ndper = 4 - endif - - ifdipolel = 0 - ifchargel = 0 - - if (ifstoklet .eq. 1) ifchargel = 1 - if (ifstrslet .eq. 1) ifdipolel = 1 - - ndper = 4 - ndl = ndper*nd - - ifpghl = 3 - ifpghtargl = 3 - -c allocate necessary arrays - - allocate(charge(ndper,nd,nsource),dipvec(ndper,nd,3,nsource), - 1 potl(ndper,nd,nsource),pottargl(ndper,nd,ntarg), - 2 gradl(ndper,nd,3,nsource),gradtargl(ndper,nd,3,ntarg), - 3 hessl(ndper,nd,6,nsource),hesstargl(ndper,nd,6,ntarg), - 4 stat=ier) - if(ier .ne. 0) then - print *, "In stfmm3d: cannot allocate Laplace call storage" - print *, "ndper =",ndper - print *, "nd =",nd - print *, "nsource =",nsource - print *, "ntarg =",ntarg - stop - endif - -c set-up appropriate vector charge and dipole arrays - -c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,l,sigma,dmu,dnu) -c$OMP$ PRIVATE(pl,pv) - do i = 1,nsource - - do j = 1,nd - do l = 1,ndper - charge(l,j,i) = 0 - dipvec(l,j,1,i) = 0 - dipvec(l,j,2,i) = 0 - dipvec(l,j,3,i) = 0 - enddo - enddo - - do j = 1,nd - if (ifstoklet .eq. 1) then - sigma(1) = stoklet(j,1,i) - sigma(2) = stoklet(j,2,i) - sigma(3) = stoklet(j,3,i) - endif - if (ifstrslet .ge. 1) then - dmu(1) = strslet(j,1,i) - dmu(2) = strslet(j,2,i) - dmu(3) = strslet(j,3,i) - dnu(1) = strsvec(j,1,i) - dnu(2) = strsvec(j,2,i) - dnu(3) = strsvec(j,3,i) - endif - - do l = 1,3 - - if (ifstoklet .eq. 1) then - charge(l,j,i) = charge(l,j,i) + sigma(l)/2 - endif - if (ifstrslet .eq. 1) then - dipvec(l,j,1,i) = dipvec(l,j,1,i) + (dmu(l)*dnu(1) + - 1 dmu(1)*dnu(l))/2 - dipvec(l,j,2,i) = dipvec(l,j,2,i) + (dmu(l)*dnu(2) + - 1 dmu(2)*dnu(l))/2 - dipvec(l,j,3,i) = dipvec(l,j,3,i) + (dmu(l)*dnu(3) + - 1 dmu(3)*dnu(l))/2 - endif - enddo - - l = 4 - - if (ifstoklet .eq. 1) then - pl = sigma(1)*source(1,i) + sigma(2)*source(2,i) + - 1 sigma(3)*source(3,i) - charge(l,j,i) = charge(l,j,i) + pl/2 - endif - if (ifstrslet .eq. 1) then - pl = dmu(1)*source(1,i) + dmu(2)*source(2,i) + - 1 dmu(3)*source(3,i) - pv = dnu(1)*source(1,i) + dnu(2)*source(2,i) + - 1 dnu(3)*source(3,i) - - dipvec(l,j,1,i) = dipvec(l,j,1,i) + - 1 (dmu(1)*pv + dnu(1)*pl)/2 - dipvec(l,j,2,i) = dipvec(l,j,2,i) + - 1 (dmu(2)*pv + dnu(2)*pl)/2 - dipvec(l,j,3,i) = dipvec(l,j,3,i) + - 1 (dmu(3)*pv + dnu(3)*pl)/2 - endif - - enddo - - enddo -c$OMP END PARALLEL DO - - -c call Laplace FMM - - iper = 0 - ier = 0 - call lfmm3d(ndl,eps,nsource,source,ifchargel,charge, - 1 ifdipolel,dipvec,iper,ifpghl,potl,gradl,hessl,ntarg, - 2 targ,ifpghtargl,pottargl,gradtargl,hesstargl,ier) - - - npt = ntarg + nsource - -c unpack stacked Laplace FMM calls - -c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,l,ifppreg1,ii) -c$OMP$ PRIVATE(pt,pl,gl,hl,vel,velgrad,press) - do i = 1,npt - - do j = 1,nd - - vel(1) = 0 - vel(2) = 0 - vel(3) = 0 - velgrad(1,1) = 0 - velgrad(2,1) = 0 - velgrad(3,1) = 0 - velgrad(1,2) = 0 - velgrad(2,2) = 0 - velgrad(3,2) = 0 - velgrad(1,3) = 0 - velgrad(2,3) = 0 - velgrad(3,3) = 0 - press = 0 - - do l = 1,ndper - - if (i .gt. ntarg) then - ifppreg1 = ifppreg - ii = i-ntarg - pt(1) = source(1,ii) - pt(2) = source(2,ii) - pt(3) = source(3,ii) - pl = potl(l,j,ii) - gl(1) = gradl(l,j,1,ii) - gl(2) = gradl(l,j,2,ii) - gl(3) = gradl(l,j,3,ii) - hl(1) = hessl(l,j,1,ii) - hl(2) = hessl(l,j,2,ii) - hl(3) = hessl(l,j,3,ii) - hl(4) = hessl(l,j,4,ii) - hl(5) = hessl(l,j,5,ii) - hl(6) = hessl(l,j,6,ii) - else - ifppreg1 = ifppregtarg - ii = i - pt(1) = targ(1,ii) - pt(2) = targ(2,ii) - pt(3) = targ(3,ii) - pl = pottargl(l,j,ii) - gl(1) = gradtargl(l,j,1,ii) - gl(2) = gradtargl(l,j,2,ii) - gl(3) = gradtargl(l,j,3,ii) - hl(1) = hesstargl(l,j,1,ii) - hl(2) = hesstargl(l,j,2,ii) - hl(3) = hesstargl(l,j,3,ii) - hl(4) = hesstargl(l,j,4,ii) - hl(5) = hesstargl(l,j,5,ii) - hl(6) = hesstargl(l,j,6,ii) - endif +c D_{ijk}(x,y) = -\delta_{jk} r_i/(4\pi r^3) - \delta_{ik} r_j/(4\\pi r^3) + \delta_{ij} r_k/(4\pi r^3) + 3 r_i r_j r_k/ (4\pi r^5) +c L_{jk} = -1/(2\pi) \delta_{jk} + 3 r_j r_k/(2\pi r^5) - - if (l .ge. 1 .and. l .le. 3) then - - vel(l) = vel(l) + pl - vel(1) = vel(1) - pt(l)*gl(1) - vel(2) = vel(2) - pt(l)*gl(2) - vel(3) = vel(3) - pt(l)*gl(3) - press = press - gl(l)*2 - if (ifppreg1 .eq. 3) then - velgrad(1,l) = velgrad(1,l) + gl(1) - velgrad(2,l) = velgrad(2,l) + gl(2) - velgrad(3,l) = velgrad(3,l) + gl(3) - - velgrad(l,1) = velgrad(l,1) - gl(1) - velgrad(l,2) = velgrad(l,2) - gl(2) - velgrad(l,3) = velgrad(l,3) - gl(3) - -c confirm hessian ordering convention... - velgrad(1,1) = velgrad(1,1) - pt(l)*hl(1) - velgrad(2,1) = velgrad(2,1) - pt(l)*hl(4) - velgrad(3,1) = velgrad(3,1) - pt(l)*hl(5) - velgrad(1,2) = velgrad(1,2) - pt(l)*hl(4) - velgrad(2,2) = velgrad(2,2) - pt(l)*hl(2) - velgrad(3,2) = velgrad(3,2) - pt(l)*hl(6) - velgrad(1,3) = velgrad(1,3) - pt(l)*hl(5) - velgrad(2,3) = velgrad(2,3) - pt(l)*hl(6) - velgrad(3,3) = velgrad(3,3) - pt(l)*hl(3) - endif - - else if (l .eq. 4) then - - vel(1) = vel(1) + gl(1) - vel(2) = vel(2) + gl(2) - vel(3) = vel(3) + gl(3) - if (ifppreg1 .eq. 3) then -c confirm hessian ordering convention... - velgrad(1,1) = velgrad(1,1) + hl(1) - velgrad(2,1) = velgrad(2,1) + hl(4) - velgrad(3,1) = velgrad(3,1) + hl(5) - velgrad(1,2) = velgrad(1,2) + hl(4) - velgrad(2,2) = velgrad(2,2) + hl(2) - velgrad(3,2) = velgrad(3,2) + hl(6) - velgrad(1,3) = velgrad(1,3) + hl(5) - velgrad(2,3) = velgrad(2,3) + hl(6) - velgrad(3,3) = velgrad(3,3) + hl(3) - endif - - endif - enddo - - if (i .gt. ntarg) then - if (ifppreg1 .ge. 1) then - pot(j,1,ii) = vel(1) - pot(j,2,ii) = vel(2) - pot(j,3,ii) = vel(3) - endif - if (ifppreg1 .ge. 2) then - pre(j,ii) = press - endif - if (ifppreg1 .ge. 3) then - grad(j,1,1,ii) = velgrad(1,1) - grad(j,2,1,ii) = velgrad(2,1) - grad(j,3,1,ii) = velgrad(3,1) - grad(j,1,2,ii) = velgrad(1,2) - grad(j,2,2,ii) = velgrad(2,2) - grad(j,3,2,ii) = velgrad(3,2) - grad(j,1,3,ii) = velgrad(1,3) - grad(j,2,3,ii) = velgrad(2,3) - grad(j,3,3,ii) = velgrad(3,3) - endif - else - if (ifppreg1 .ge. 1) then - pottarg(j,1,ii) = vel(1) - pottarg(j,2,ii) = vel(2) - pottarg(j,3,ii) = vel(3) - endif - if (ifppreg1 .ge. 2) then - pretarg(j,ii) = press - endif - if (ifppreg1 .ge. 3) then - gradtarg(j,1,1,ii) = velgrad(1,1) - gradtarg(j,2,1,ii) = velgrad(2,1) - gradtarg(j,3,1,ii) = velgrad(3,1) - gradtarg(j,1,2,ii) = velgrad(1,2) - gradtarg(j,2,2,ii) = velgrad(2,2) - gradtarg(j,3,2,ii) = velgrad(3,2) - gradtarg(j,1,3,ii) = velgrad(1,3) - gradtarg(j,2,3,ii) = velgrad(2,3) - gradtarg(j,3,3,ii) = velgrad(3,3) - endif - endif - enddo - enddo -c$OMP END PARALLEL DO - - return - end - -c********************************************************************** -c -c We take the following conventions for the Stokes kernels. -c -c 1) The dynamic viscosity (mu) is assumed to be 1. -c 2) All kernels are using standard definitions. -c -c For a source y and target x, let r_i = x_i-y_i (note sign) -c and let r = sqrt(r_1^2 + r_2^2 + r_3^2) -c -c The Stokeslet, G_{ij}, and its associated pressure tensor, P_j, -c we define as -c -c G_{ij}(x,y) = ( delta_{ij}/r + r_i r_j / r^3 )/(8\pi) -c P_j(x,y) = 1/(4\pi) * r_j/r^3 -c -c The (Type I) stresslet, T_{ijk}, and its associated pressure -c tensor, PI_{jk}, we define as -c -c T_{ijk}(x,y) = 3/(4\pi) r_i r_j r_k/ r^5 -c PI_{jk}(x,y) = -1/(2\pi) delta_{jk}/r^3 + 3/(2\pi) r_j r_k/r^5 - subroutine stfmm3d_new(nd, eps, + subroutine stfmm3d(nd, eps, $ nsource, source, $ ifstoklet, stoklet, ifstrslet, strslet, strsvec, $ ifrotlet, rotlet, rotvec, @@ -514,10 +62,14 @@ subroutine stfmm3d_new(nd, eps, c c u_i(x) = sum_m G_{ij}(x,y^{(m)}) sigma^{(m)}_j c + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k +c + sum_m R_{ijk}(x,y^{(m)}) rlet^{(m)}_j rvec^{(m)}_k +c + sum_m D_{ijk}(x,y^{(m)}) dlet^{(m)}_j dvec^{(m)}_k c c for x at all of the target locations, and i=1,2,3, c where sigma^{(m)} is the Stokeslet force, mu^{(m)} is the -c stresslet strength, and nu^{(m)} is the stresslet orientation +c stresslet strength, nu^{(m)} is the stresslet orientation, +c rlet^{(m)} is the rotlet strength, rvec^{(m)} is the rotlet +c dlet^{(m)} is the doublet strength, and dvec^{(m)} is the doublet c (note that each of these is a 3 vector per source point y^{(m)}). c Repeated indices are taken as summed over 1,2,3, ie, Einstein c convention. For x a source point, the self-interaction in the @@ -528,9 +80,12 @@ subroutine stfmm3d_new(nd, eps, c c p(x) = sum_m P_j(x,y^m) sigma^{(m)}_j c + sum_m PI_{jk}(x,y{(m)}) mu^{(m)}_j nu^{(m)}_k +c + sum_m L_{jk}(x,y^{(m)}) dlet^{(m)}_j dvec^{(m)}_k c c grad_l u_i(x) = grad_l [sum_m G_{ij}(x,y^m) sigma^{(m)}_j c + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k] +c + sum_m R_{ijk}(x,y^{(m)}) rlet^{(m)}_j rvec^{(m)}_k +c + sum_m D_{ijk}(x,y^{(m)}) dlet^{(m)}_j dvec^{(m)}_k c c Note that these two may be combined to get the stress tensor. c diff --git a/src/Stokes/stokkernels.f b/src/Stokes/stokkernels.f index 9844299..1c6d459 100644 --- a/src/Stokes/stokkernels.f +++ b/src/Stokes/stokkernels.f @@ -20,6 +20,16 @@ c T_{ijk}(x,y) = 3/(4\pi) r_i r_j r_k/ r^5 c PI_{jk} = -1/(2\pi) delta_{jk} + 3/(2\pi) r_j r_k/r^5 c +c The rotlet, R_{ijk}, and its associated pressure tensor, Q_{jk}, are +c +c R_{ijk}(x,y) = -\delta_{ik} r_j/(4\pi r^3) + \delta_{ij} r_k/(4\pi r^3) +c Q_{jk} = 0; +c +c The doublet, D_{ijk}, and its associated pressure tensor, L_{jk}, are +c +c D_{ijk}(x,y) = -\delta_{jk} r_i/(4\pi r^3) - \delta_{ik} r_j/(4\\pi r^3) + \delta_{ij} r_k/(4\pi r^3) + 3 r_i r_j r_k/ (4\pi r^5) +c L_{jk} = -1/(2\pi) \delta_{jk} + 3 r_j r_k/(2\pi r^5) + subroutine st3ddirectstokg(nd,sources,stoklet,ns,targ,nt, 1 pot,pre,grad,thresh) @@ -423,15 +433,22 @@ subroutine st3ddirectstokstrsrotdoubg(nd,sources,stoklet,istress, c c pot(x) = pot(x) + sum_m G_{ij}(x,y^{(m)}) sigma^{(m)}_j c + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k +c + sum_m R_{ijk}(x,y^{(m)}) rlet^{(m)}_j rvec^{(m)}_k +c + sum_m D_{ijk}(x,y^{(m)}) dlet^{(m)}_j dvec^{(m)}_k c c pre(x) = sum_m P_j(x,y^m) sigma^{(m)}_j -c + sum_m T_{ijk}(x,y^{(m)}) PI_{jk} mu^{(m)}_j nu^{(m)}_k +c + sum_m PI_{jk} mu^{(m)}_j nu^{(m)}_k +c + sum_m L_{jk}(x,y^{(m)}) dlet^{(m)}_j dvec^{(m)}_k c c grad(x) = Grad[sum_m G_{ij}(x,y^m) sigma^{(m)}_j c + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k] +c + sum_m R_{ijk}(x,y^{(m)}) rlet^{(m)}_j rvec^{(m)}_k +c + sum_m D_{ijk}(x,y^{(m)}) dlet^{(m)}_j dvec^{(m)}_k c c where sigma^{(m)} is the Stokeslet charge, mu^{(m)} is the -c stresslet charge, and nu^{(m)} is the stresslet orientation +c stresslet charge, nu^{(m)} is the stresslet orientation +c rlet^{(m)} is the rotlet strength, rvec^{(m)} is the rotlet +c dlet^{(m)} is the doublet strength, and dvec^{(m)} is the doublet c (note that each of these is a 3 vector per source point y^{(m)}). c For x a source point, the self-interaction in the sum is omitted. c diff --git a/test/Stokes/test_stfmm3d.f b/test/Stokes/test_stfmm3d.f index 7462894..0f82b9f 100644 --- a/test/Stokes/test_stfmm3d.f +++ b/test/Stokes/test_stfmm3d.f @@ -14,6 +14,8 @@ real *8, allocatable :: stoklet(:,:), strslet(:,:) real *8, allocatable :: strsvec(:,:) + real *8, allocatable :: rotlet(:,:), rotvec(:,:) + real *8, allocatable :: doublet(:,:), doubvec(:,:) integer *8 ipass(10) integer *8 nd,ns,nt,ifstoklet,ifstrslet,ifppreg,ifppregtarg,ier @@ -70,14 +72,17 @@ eps = 1d-9 ifstoklet = 1 ifstrslet = 1 + ifrotlet = 0 + ifdoublet = 0 ifppreg = 3 ifppregtarg = 3 call cpu_time(t1) c$ t1 = omp_get_wtime() call stfmm3d(nd,eps,ns,source,ifstoklet,stoklet, - 1 ifstrslet,strslet,strsvec,ifppreg,pot,pre,grad, - 2 nt,targ,ifppregtarg,pottarg,pretarg,gradtarg,ier) + 1 ifstrslet,strslet,strsvec,ifrotlet,rotlet,rotvec, + 2 ifdoublet,doublet,doubvec,ifppreg,pot,pre,grad, + 3 nt,targ,ifppregtarg,pottarg,pretarg,gradtarg,ier) call cpu_time(t2) c$ t2 = omp_get_wtime() diff --git a/test/Stokes/test_stokkernels_rotlet_doublet.f b/test/Stokes/test_stokkernels_rotlet_doublet.f index bb000ee..b73ae8b 100644 --- a/test/Stokes/test_stokkernels_rotlet_doublet.f +++ b/test/Stokes/test_stokkernels_rotlet_doublet.f @@ -234,13 +234,13 @@ errmax = 0 ifppgreg = 1 ifppregtarg = 3 - call stfmm3d_new(nd1,1d-12,ns,source,istoklet,stoklet, + call stfmm3d(nd1,1d-12,ns,source,istoklet,stoklet, 1 istress,strslet,strsvec, 2 irotlet,rotstr,rotvec, 3 idoublet,doubstr,doubvec, 4 ifppgreg,pot_src,pre_src,grad_src, 5 nt,targ,ifppregtarg,pot3,pre3,grad3,ier) - call stfmm3d_new(nd1,1d-12,ns,source,istoklet,stoklet, + call stfmm3d(nd1,1d-12,ns,source,istoklet,stoklet, 1 istress,strslet,strsvec, 2 irotlet,rotstr,rotvec, 3 idoublet,doubstr,doubvec,