diff --git a/docs/math.rst b/docs/math.rst index 294b683..d339bda 100644 --- a/docs/math.rst +++ b/docs/math.rst @@ -1,7 +1,7 @@ Definitions =========== -Let $x_{j} \in \mathbb{R}^{2}$, $i=1,2,\ldots N$, denote a collection -of source locations and let $t_{i} \in \mathbb{R}^{2}$ denote a collection +Let $x_{j} = (x_{j,1}, x_{j,2}) \in \mathbb{R}^{2}$, $i=1,2,\ldots N$, denote a collection +of source locations and let $t_{i} = (t_{i,1}, t_{i,2}) \in \mathbb{R}^{2}$ denote a collection of target locations. @@ -13,12 +13,12 @@ The Laplace FMM comes in three varieties gradients, and hessians are all double precision * lfmm2d: All of the above quantities are double complex * cfmm2d: Same as lfmm2d except that there are no dipole orientation vectors, - and the dipole kernel is given by $1/(zt_{i} - zx_{j})$ where $zt_{i}, zx_{j}$ + and the dipole kernel is given by $1/(z_{i} - \xi_{j})$ where $z_{i}, \xi_{j}$ are the complex versions of the source and target locations given by - $zt_{i} = t_{i}(1) + i*t_{i}(2)$, and $zx_{j} = x_{j}(1) + i*y_{j}(2)$. + $z_{i} = t_{i,1} + i\cdot t_{i,2}$, and $\xi_{j} = x_{j,1} + i \cdot x_{j,2}$. Moreover, the gradients, and hessians in this case are derivatives with respect to $z$, and in a slight abuse of - notation, we set $d/dz log(|z|) = 1/z$ + notation, we set $\frac{\mathrm{d}}{\mathrm{d} z} \log(\|z\|) = 1/z$ Laplace FMM (rfmm2d) @@ -31,14 +31,14 @@ denote a collection of dipole strengths, and $d_{j} \in \mathbb{R}^{2}$, $j=1,2,\ldots N$, denote the corresponding dipole orientation vectors. The Laplace FMM (rfmm2d) computes -the potential $u(x) \in \mathbb{R}$ and the its gradient $\nabla u(x) \in +the potential $u(x) \in \mathbb{R}$ and its gradient $\nabla u(x) \in \mathbb{R}^{2}$ given by .. math:: :label: rlap_nbody - u(x) = \sum_{j=1}^{N} c_{j} log(\|x-x_{j}\|) - v_{j} d_{j}\cdot \nabla log(\|x-x_{j}\|) \, , + u(x) = \sum_{j=1}^{N} c_{j} \log{(\|x-x_{j}\|)} - v_{j} d_{j}\cdot \nabla \log{(\|x-x_{j}\|)} \, , at the source and target locations. When $x=x_{j}$, the term corresponding to $x_{j}$ is dropped from the sum. @@ -54,14 +54,14 @@ denote a collection of dipole strengths, and $d_{j} \in \mathbb{R}^{2}$, $j=1,2,\ldots N$, denote the corresponding dipole orientation vectors. The Laplace FMM (rfmm2d) computes -the potential $u(x) \in \mathbb{C}$ and the its gradient $\nabla u(x) \in +the potential $u(x) \in \mathbb{C}$ and its gradient $\nabla u(x) \in \mathbb{C}^{2}$ given by .. math:: :label: llap_nbody - u(x) = \sum_{j=1}^{N} c_{j} log(\|x-x_{j}\|) - v_{j} d_{j}\cdot \nabla log(\|x-x_{j}\|) \, , + u(x) = \sum_{j=1}^{N} c_{j} \log{(\|x-x_{j}\|)} - v_{j} d_{j}\cdot \nabla \log{(\|x-x_{j}\|)} \, , at the source and target locations. When $x=x_{j}$, the term corresponding to $x_{j}$ is dropped from the sum. @@ -76,18 +76,18 @@ $j=1,2,\ldots N$, denote a collection of dipole strengths. The Laplace FMM (cfmm2d) computes -the potential $u(z) \in \mathbb{C}$ and the its gradient $d/dz u(z) \in +the potential $u(z) \in \mathbb{C}$ and its gradient $\frac{\mathrm{d}}{\mathrm{d}z} u(z) \in \mathbb{C}$ given by .. math:: :label: clap_nbody - u(z) = \sum_{j=1}^{N} c_{j} log(\|z-zx_{j}\|) - \frac{v_{j}}{z-zx_{j}} \, , + u(z) = \sum_{j=1}^{N} c_{j} \log{(\|z-\xi_{j}\|)} - \frac{v_{j}}{z-\xi_{j}} \, , -at the source and target locations. When $z=zx_{j}$, the term -corresponding to $zx_{j}$ is dropped from the sum. Here -$z = x(1) + i x(2)$, and $zx_{j} = x_{j}(1) + ix_{j}(2)$, are the complex +at the source and target locations. When $z=\xi_{j}$, the term +corresponding to $\xi_{j}$ is dropped from the sum. Here +$z = x(1) + i x(2)$, and $\xi_{j} = x_{j}(1) + ix_{j}(2)$, are the complex numbers corresponding to $x$ and $x_{j}$ respectively. @@ -103,7 +103,7 @@ Let $k\in\mathbb{C}$ denote the wave number or the Helmholtz parameter. The Helmholtz FMM computes -the potential $u(x)$ and the its gradient $\nabla u(x)$ +the potential $u(x)$ and its gradient $\nabla u(x)$ given by .. math:: @@ -118,13 +118,155 @@ corresponding to $x_{j}$ is dropped from the sum. Biharmonic FMM *************** +Let $c_{j} = (c_{j,1}, c_{j,2})\in \mathbb{C}^2$, +$j=1,2,\ldots N$, +denote a collection of charge strengths, and +$v_{j} = (v_{j,1}, v_{j,2}, v_{j,3}) \in \mathbb{C}^{3}$, +$j=1,2,\ldots N$, +denote a collection of dipole strengths. + +The Biharmonic FMM computes +the potential $u(x)$ and its `gradient` = $(P_{z} \frac{\mathrm{d}}{\mathrm{d}z}, P_{\overline{z}} \frac{\mathrm{d}}{\mathrm{d}z}, \frac{\mathrm{d}}{\mathrm{d}\overline{z}})$ +given by + +.. math:: + :label: biharm_nbody + + u(z) &= \sum_{j=1}^{N} c_{j,1} \log{\|z - \xi_{j}\|} + + c_{j,2} \frac{z-\xi_{j}}{\overline{z-\xi_{j}}} + \frac{v_{j,1}}{z - \xi_{j}} + + \frac{v_{j,3}}{\overline{z-\xi_{j}}} + + v_{j,2} \frac{z - \xi_{j}}{(\overline{z-\xi_{j}})^2} \, , \\ + P_{z} \frac{\mathrm{d}}{\mathrm{d} z}u(z) &= \sum_{j=1}^{N} \frac{c_{j,1}}{z-\xi_{j}} - + \frac{v_{j,1}}{(z-\xi_{j})^2} \, \\ + P_{\overline{z}} \frac{\mathrm{d}}{\mathrm{d} z} u(z) &= + \sum_{j=1}^{N} \frac{c_{j,2}}{\overline{z-\xi_{j}}} + + \frac{v_{j,2}}{(\overline{z-\xi_{j}})^2} \, + ,\\ + \frac{\mathrm{d}}{\mathrm{d}\overline{z}} u(z) &= + \sum_{j=1}^{N} \frac{c_{j,1}}{\overline{z-\xi_{j}}} - + c_{j,2} \frac{z-\xi_{j}}{(\overline{z-\xi_{j}})^2} - + \frac{v_{j,3}}{(\overline{z-\xi_{j}})^2} - + 2v_{j,2} \frac{z - \xi_{j}}{(\overline{z-\xi_{j}})^3} \, , \\ + +at the source and target locations. When $z=\xi_{j}$, the term +corresponding to $\xi_{j}$ is dropped from the sum. Here +$z = x(1) + i x(2)$, and $\xi_{j} = x_{j}(1) + ix_{j}(2)$, are the complex +numbers corresponding to $x$ and $x_{j}$ respectively. +When $x=x_{j}$, the term +corresponding to $x_{j}$ is dropped from the sum. + Modified Biharmonic FMM ************************ +Let $G^{\textrm{mbh}}(x,y)$ denote the modified biharmonic +Green's function given by + +.. math:: + G^{\textrm{mbh}}(x,y) = \frac{1}{2\pi \beta^2}\left(K_{0}(\beta \|x-y \|) - \log{\|x-y\|}\right) + +where $K_{0}$ is the modified Bessel function of order $0$, and $\beta$ is the +Modified Biharmonic wavenumber. + + +Let $c_{j} \in \mathbb{C}$, +denote a collection of charge strengths, +$v_{j} \in \mathbb{C}$, +denote a collection of dipole strengths, +$d_{j} = (d_{j,1}, d_{j,2})$ denote a collection +of dipole vectors, +$q_{j} \in \mathbb{C}$ denote a collection of +quadrupole strengths, +$w_{j} = (w_{j,1}, w_{j,2}, w_{j,3}) \in \mathbb{C}^{3}$, +denote a collection of quadrupole orientation vectors, +$o_{j}$ denote a collection of octopole strengths, and +$p_{j} = (p_{j,1}, p_{j,2}, p_{j,3}, p_{j,4}) \in \mathbb{C}^{4}$, +denote a collceetion of octopole strengths. For all the vectors, +$j=1,2,\ldots N$. + +The Modified Biharmonic FMM computes +the potential $u(x)$ and its gradients, and hessians, +given by + +.. math:: + :label: modbiharm_nbody + + u(x) = \sum_{j=1}^{N} &c_{j} G^{\textrm{mbh}}(x,x_{j}) + + v_{j} d_{j} \cdot \nabla_{y} G^{\textrm{mbh}}(x,x_{j}) + \\ + &q_{j} \left(w_{j,1} \partial_{y_{1},y_{1}} G^{\textrm{mbh}}(x,x_{j}) + w_{j,2} + \partial_{y_{1},y_{2}} G^{\textrm{mbh}}(x,x_{j}) + w_{j,3} + \partial_{y_{2},y_{2}} G^{\textrm{mbh}}(x,x_{j}) \right) + \\ + &o_{j} \big( p_{j,1} \partial_{y_{1},y_{1},y_{1}} G^{\textrm{mbh}}(x,x_{j}) + + p_{j,2} \partial_{y_{1},y_{1},y_{2}} G^{\textrm{mbh}}(x,x_{j}) + \\ + &p_{j,3} \partial_{y_{1},y_{2},y_{2}} G^{\textrm{mbh}}(x,x_{j}) + + p_{j,4} \partial_{y_{2},y_{2},y_{2}} G^{\textrm{mbh}}(x,x_{j}) \big) \, , + +at the source and target locations. When $x=x_{j}$, the term +corresponding to $x_{j}$ is dropped from the sum. + Stokes FMM ************ +Let $G^{\textrm{stok}}(x,y)$ denote the Stokeslet given by + +.. math:: + G^{\textrm{stok}}(x,y) = \frac{1}{2} + \begin{bmatrix} + -\log{\|x-y \|} + \frac{(x_{1}-y_{1})^2}{\|x-y\|^2} & \frac{(x_{1}-y_{1}) + (x_{2}-y_{2})}{\|x-y \|^2} \\ + \frac{(x_{1}-y_{1})(x_{2}-y_{2})}{\|x-y \|^2} & + -\log{\|x-y \|} + \frac{(x_{2}-y_{2})^2}{\|x-y \|^2} + \end{bmatrix} \, , + +$P^{\textrm{stok}}(x,y)$ denote the associated pressure tensor + +.. math:: + P(x,y) = \frac{1}{\|x-y \|^2}\begin{bmatrix} + (x_{1}-y_{1}) & + (x_{2}-y_{2}) + \end{bmatrix} \, , + +$T^{\textrm{stok}}(x,y)$ denote the Stresslet whose action of a vector $v$ +is given by + +.. math:: + v\cdot T^{\textrm{stok}}(x,y) = -\frac{2 v \cdot (x-y)}{\|x-y \|^4} + \begin{bmatrix} + (x_{1} - y_{1})^2 & (x_{1}-y_{1})(x_{2}-y_{2}) \\ + (x_{1}-y_{1})(x_{2}-y_{2}) & (x_{2} - y_{2})^2 + \end{bmatrix} \, , + +and finally let $\Pi^{\textrm{stok}} (x,y)$ denote its assocaited pressure tensor given by + +.. math:: + v\cdot \Pi(x,y)^{textrm{stok}} = -\frac{v}{\|x-y \|^2} + \frac{2 v \cdot(x-y)}{\|x-y \|^4} + \begin{bmatrix} + (x_{1}-y_{1}) \\ + (x_{2}-y_{2}) + \end{bmatrix} + +Let $c_{j} \in \mathbb{R}^2$, +$j=1,2,\ldots N$, +denote a collection of stokeslet strengths, $v_{j} \in \mathbb{R}^2$, +$j=1,2,\ldots N$, +denote a collection of stresslet strengths, and $d_{j} \in \mathbb{R}^{2}$, +$j=1,2,\ldots N$, denote the corresponding stresslet orientation vectors. + +The Stokes FMM computes +the potential $u(x) \in \mathbb{R}^2$, its gradient $\nabla u(x) \in +\mathbb{R}^{2\times 2}$, and the pressure $p$ given by + +.. math:: + :label: stok_nbody + + u(x) &= \sum_{j=1}^{N} G^{\textrm{stok}}(x,x_{j}) c_{j} + d_{j} \cdot + T^{\textrm{stok}}(x,x_{j}) \cdot v_{j} \, , \\ + p(x) &= \sum_{j=1}^{N} c_{j} P^{\textrm{stok}}(x,x_{j}) + d_{j} \cdot + \Pi^{\textrm{stok}}(x,x_{j}) \cdot v_{j}^{T} + +at the source and target locations. When $x=x_{j}$, the term +corresponding to $x_{j}$ is dropped from the sum. + Vectorized versions ******************* The vectorized versions of the Helmholtz FMM, diff --git a/python/bhfmmexample.py b/python/bhfmmexample.py new file mode 100644 index 0000000..ae25597 --- /dev/null +++ b/python/bhfmmexample.py @@ -0,0 +1,34 @@ +#!/usr/bin/env python + +import fmm2dpy as fmm +import numpy as np + + +# +# This is a sample code to demonstrate how to use +# the fmm libraries +# + +# sample with one density, sources to sources, +# charge interactions, and potential only +# +n = 100000 +nd = 1 +sources = np.random.uniform(0,1,(2,n)) +eps = 10**(-5) + +charges = np.random.uniform(0,1,(2,n)) + 1j*np.random.uniform(0,1,(2,n)) +dipoles = np.random.uniform(0,1,(3,n)) + 1j*np.random.uniform(0,1,(3,n)) +out = fmm.bhfmm2d(eps=eps,sources=sources,charges=charges,dipoles=dipoles,pg=1) + + +# sample with a vector of densities, sources to +# sources and targets, dipole interactions, +# potential and gradietns + +nd = 3 +nt = 1870 +targ = np.random.uniform(0,1,(2,nt)) +dipoles = np.random.uniform(0,1,(nd,3,n)) + 1j*np.random.uniform(0,1,(nd,3,n)) +out2 = fmm.bhfmm2d(eps=eps,sources=sources,dipoles=dipoles,\ + targets=targ,nd=nd,pg=2,pgt=2) diff --git a/src/biharmonic/bhfmm2d.f b/src/biharmonic/bhfmm2d.f index a552451..6182a0d 100644 --- a/src/biharmonic/bhfmm2d.f +++ b/src/biharmonic/bhfmm2d.f @@ -51,7 +51,7 @@ subroutine bhfmm2d(nd,eps,ns,sources,ifcharge,charge, c and velocity are interchangably used c c -c vel(z) = \sum charge1_j *\log|z-z_j| + +c vel(z) = \sum 2*charge1_j *\log|z-z_j| + c (charge2_j) (z-z_j)/(z-z_j)_bar+ c dip1_j/(z-z_j) + dip2_j/(z-z_j)_bar- c dip3_j (z-z_j)/(z-z_j)^2_bar @@ -65,8 +65,8 @@ subroutine bhfmm2d(nd,eps,ns,sources,ifcharge,charge, c \phi(z) = \sum charge_j log(z-z_j)+ c dip1/(z-z_j) c -c \psi(z) = \sum dippar2_bar/(z-z_j)+ -c z_j_bar dippar1/(z-z_j)^2 + +c \psi(z) = \sum dip2_bar/(z-z_j)+ +c z_j_bar dip1/(z-z_j)^2 + c c In all of the sums above, the terms for which |z-z_{j}| <= 2^-51|b| c where |b| is the size of the root box will be ignored. diff --git a/src/modified-biharmonic/mbhfmm2d.f b/src/modified-biharmonic/mbhfmm2d.f index 57cb858..3c3f786 100644 --- a/src/modified-biharmonic/mbhfmm2d.f +++ b/src/modified-biharmonic/mbhfmm2d.f @@ -36,10 +36,10 @@ subroutine mbhfmm2d(nd,eps,beta,ns,sources,ifcharge,charge, c u(x) = sum_j charge(j)*G(x,y(j)) c + dipstr(j)*(G_{y1}(x,y(j))*dipvec(1,j) + G_{y2}*dipvec(2,j)) c + quadstr(j)*(G_{y1y1}(x,y(j))*quadvec(1,j) -c + G_{y1y2}*quadvec(2,j) + G_{y2y2}*quadvec(2,j)) +c + G_{y1y2}*quadvec(2,j) + G_{y2y2}*quadvec(3,j)) c + octstr(j)*(G_{y1y1y1}(x,y(j))*octvec(1,j) -c + G_{y1y1y2}*octvec(2,j) + G_{y1y2y2}*octvec(2,j) -c + G_{y1y2y2}*octvec(2,j)) +c + G_{y1y1y2}*octvec(2,j) + G_{y1y2y2}*octvec(3,j) +c + G_{y2y2y2}*octvec(4,j)) c c INPUT PARAMETERS: c diff --git a/src/stokes/stfmm2d.f b/src/stokes/stfmm2d.f index 744e9bc..f8f9d3a 100644 --- a/src/stokes/stfmm2d.f +++ b/src/stokes/stfmm2d.f @@ -11,16 +11,16 @@ c The Stokeslet, G_{ij}, and its associated pressure tensor, P_j, c (without the 1/2pi scaling) are c -c UPDATE formulae -c -c G_{ij}(x,y) = (r_i r_j)/(2r^3) + delta_{ij}/(2r) -c P_j(x,y) = r_j/r^3 +c G_{ij}(x,y) = (r_i r_j)/(2r^2) - delta_{ij}log(r)/(2) +c P_j(x,y) = r_j/r^2 c c The (Type I) stresslet, T_{ijk}, and its associated pressure -c tensor, PI_{jk}, (without the 1/4pi scaling) are +c tensor, PI_{jk}, (without the 1/2pi scaling) are c -c T_{ijk}(x,y) = -3 r_i r_j r_k/ r^5 -c PI_{jk} = -2 delta_{jk} + 6 r_j r_k/r^5 +c T_{ijk}(x,y) = -2 r_i r_j r_k/ r^4 +c PI_{jk} = - delta_{jk}/r^2 + 2 r_j r_k/r^4 +c +c c