Merge pull request #9 from mipals/adding_laplace_wrappers

Adding Laplace wrappers

.github/workflows/CI.yml

@@ -18,7 +18,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.8'
+          - '1.6'
+          - '1.10'
           - 'nightly'
         os:
           - ubuntu-latest
@@ -34,3 +35,7 @@ jobs:
       - uses: julia-actions/cache@v1
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v4
+        env:
+          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}

.gitignore

@@ -1,3 +1,4 @@
 
 .DS_Store
 Manifest.toml
+test/tests.jl

Project.toml

@@ -1,22 +1,28 @@
 name = "FMM2D"
 uuid = "2d63477d-9690-4b75-bcc1-c3461d43fecc"
 authors = ["Mikkel Paltorp Schmitt"]
-version = "0.1.0"
+version = "0.2.0"
 
 [deps]
 FMM2D_jll = "0fc7e017-fbf9-5352-9b8d-9e37f15ce1cd"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
+FMM2D_jll = "1.1"
+Aqua = "0.8"
 SafeTestsets = "0.1.0"
 SpecialFunctions = "1.8.1,2"
-julia = "1.6,1.7,1.8,1.9"
+LinearAlgebra = "1"
+Test = "1"
+Random = "1"
+julia = "1"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [targets]
-test = ["Test"]
+test = ["Test", "Aqua", "LinearAlgebra", "Random", "SafeTestsets", "SpecialFunctions"]

README.md

@@ -1,10 +1,11 @@
 # FMM2D.jl
 
 [![Build Status](https://github.com/mipals/FMM2D.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/mipals/FMM2D.jl/actions/workflows/CI.yml?query=branch%3Amain)
+[![Coverage](https://codecov.io/gh/mipals/FMM2D.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/mipals/FMM2D.jl)
 
 FMM2D.jl is a Julia interface for computing N-body interactions using the [Flatiron Institute's FMM2D library](https://github.com/flatironinstitute/fmm2d/).
 
-Currently, the interface only supports Helmholtz problems.
+Currently, the wrapper only wraps the Helmholtz and Laplace functionalities.
 
 ## Helmholtz FMM
 Let $c_j \in \mathbb{C},\ j = 1,\dots,N$ denote a collection of charge strengths, $v_j \in \mathbb{C},\ j = 1,\dots,N$ denote a collection of dipole strengths, and $\mathbf{d}_j\in\mathbb{R}^2,\ j = 1,\dots,N$ denote the corresponding dipole orientation vectors. Furthermore, $k \in \mathbb{C}$ denote the wave number.
@@ -18,21 +19,116 @@ $$
 where $H_0^{(1)}$ is the Hankel function of the first kind of order 0. When $\mathbf{x} = \mathbf{x}_j$ the $j$ th term is dropped from the sum.
 
 
-## Example
+### Example
 ```julia
 using FMM2D
 
 # Simple example for the FMM2D Library
 thresh = 10.0^(-5)          # Tolerance
 zk     = rand(ComplexF64)   # Wavenumber
 
-# Source-to-source,
+# Source-to-source
 n = 200
 sources = rand(2,n)
 charges = rand(ComplexF64,n)
-vals = hfmm2d(thresh,zk,sources,charges=charges,pg=1)
+pg = 3 # Evaluate potential, gradient, and Hessian at the sources
+vals = hfmm2d(eps=thresh,zk=zk,sources=sources,charges=charges,pg=pg)
+vals.pot
+vals.grad
+vals.hess
+
+# Source-to-target
+m = 200
+targets = rand(2,m)
+pgt = 3 # Evaluate potential, gradient, and Hessian at the targets
+vals = hfmm2d(targets=targets,eps=thresh,zk=zk,sources=sources,charges=charges,pgt=pgt)
+vals.pottarg
+vals.gradtarg
+vals.hesstarg
+```
+
+## Laplace
+The Laplace problem in 2D have the following form
+
+$$
+    u(x) = \sum_{j=1}^{N} \left[c_{j} \text{log}\left(\|x-x_{j}\|\right) - d_{j}v_{j} \cdot \nabla( \text{log}(\|x-x_{j}\|) )\right],
+$$
+
+In the case of complex charges and dipole strengths ($c_j, v_j \in \mathbb{C}^n$) the function call `lfmm2d` has to be used. In the case of real charges and dipole strengths ($c_j, v_j \in \mathbb{R}^n$) the function call `rfmm2d` has to be used.
+
+
+### Example
+```julia
+using FMM2D
+
+# Simple example for the FMM2D Library
+thresh = 10.0^(-5)          # Tolerance
+
+# Source-to-source
+n = 200
+sources = rand(2,n)
+charges = rand(ComplexF64,n)
+dipvecs = randn(2,n)
+dipstr = rand(ComplexF64,n)
+pg = 3 # Evaluate potential, gradient, and Hessian at the sources
+vals = lfmm2d(eps=thresh,sources=sources,charges=charges,dipvecs=dipvecs,dipstr=dipstr,pg=pg)
+vals.pot
+vals.grad
+vals.hess
+
+# Source-to-target
+m = 100
+targets = rand(2,m)
+pgt = 3 # Evaluate potential, gradient, and Hessian at the targets
+vals = lfmm2d(targets=targets,eps=thresh,sources=sources,charges=charges,dipvecs=dipvecs,dipstr=dipstr,pgt=pgt)
+vals.pottarg
+vals.gradtarg
+vals.hesstarg
 ```
 
+
+<!-- In addition, the `FMM2D` library also includes the following sum
+
+$$
+    u(z) = \sum_{j=1}^{N} \left[c_{j} \text{log}\left(z-\varepsilon_j\right) - \frac{v_j}{z - \varepsilon_j}\right],
+$$
+
+where $c_j \ in$ -->
+
+<!-- ## Biharmonic
+
+Let $c_j = (c_{j,1}, c_{j,2}) \in \mathbb{C}^2,\ j=1,2,\dots, 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,\dots,N$ denote a collection of dipole strengths. 
+
+The Biharmonic FMM computes the potential $u(x)$ 
+
+$$
+    u(z) = \sum_{j=1}^N \left[c_{j,1}\text{log}\left(\|z - \varepsilon_j\|\right) + c_{j,2}\frac{z - \varepsilon_j}{\overline{z - \varepsilon_j}} + \frac{v_{j,1}}{z - \varepsilon_j} + \frac{v_{j,3}}{\overline{z - \varepsilon_j}} + v_{j,2}\frac{z - \varepsilon_j}{\left(\overline{z - \varepsilon_j}\right)^2}\right],
+$$
+
+as well as 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
+
+$$
+\begin{aligned}
+    P_z\frac{\mathrm{d}}{\mathrm{d}z}u(z) &= \sum_{j=1}^{N}\left[\frac{c_{j,1}}{z - \varepsilon_j} - \frac{v_{j,1}}{\left(z - \varepsilon_j\right)^2}\right]\\
+    P_{\overline{z}}\frac{\mathrm{d}}{\mathrm{d}z}u(z) &= \sum_{j=1}^{N}\left[\frac{c_{j,2}}{\overline{z - \varepsilon_j}} - \frac{v_{j,2}}{\left(\overline{z - \varepsilon_j}\right)^2}\right]\\
+    \frac{\mathrm{d}}{\mathrm{d}\overline{z}}u(z) &= \sum_{j=1}^{N}\left[\frac{c_{j,1}}{\overline{z - \varepsilon_j}} - c_{j,2}\frac{z - \varepsilon_{j}}{\left(\overline{z - \varepsilon_j}\right)^2}- v_{j,3}\frac{z - \varepsilon_{j}}{\left(\overline{z - \varepsilon_j}\right)^2} - 2v_{j,2}\frac{z - \varepsilon_{j}}{\left(\overline{z - \varepsilon_j}\right)^3}\right]
+\end{aligned}.
+$$ 
+
+at the source and target locations. When $z = \varepsilon_j$ the $j$ th term is dropped from the sum. Note here that $z = x_1 + i x_2$ and $\varepsilon_j = x_{j,1} + ix_{j,2}$. -->
+
+
+
+## Stokes
+
+$$
+    u(x) = \sum_{j=1}^NG^\text{stok}(x,x_j)c_j + d_j\cdot T^\text{stok}(x,x_j)\cdot v_j
+$$
+
+$$
+    p(x) = \sum_{j=1}^NP^\text{stok}(x,x_j)c_j + d_j\cdot \Pi^\text{stok}(x,x_j)\cdot v_j^\top
+$$
+
 ## Related Package
-[FMMLIB2D.jl](https://github.com/ludvigak/FMMLIB2D.jl) interfaces the [FMMLIB2D](https://github.com/zgimbutas/fmmlib2d) library which the [FMM2D library improves on](https://fmm2d.readthedocs.io/en/latest/). Due to missing sufficient documentation this package currently only supports Helmholtz problems whereas FMMLIB2D.jl also includes Laplace problems. A future release of this package is expected to also support Laplace problems as well as the other supported kernels from the FMM2D library.  
+[FMMLIB2D.jl](https://github.com/ludvigak/FMMLIB2D.jl) interfaces the [FMMLIB2D](https://github.com/zgimbutas/fmmlib2d) library which the [FMM2D library improves on](https://fmm2d.readthedocs.io/en/latest/). 
 

examples/hfmm2d_example.jl

@@ -8,10 +8,10 @@ zk     = rand(ComplexF64)   # Wavenumber
 n = 200
 sources = rand(2,n)
 charges = rand(ComplexF64,n)
-vals    = hfmm2d(thresh,zk,sources,charges=charges,pg=1)
+vals    = hfmm2d(eps=thresh,zk=zk,sources=sources,charges=charges,pg=1)
 
 # Source-to-target
 ntargets = 300
 targets  = rand(2,ntargets)
-vals     = hfmm2d(thresh,zk,sources;charges=charges,targets=targets,pgt=1)
+vals     = hfmm2d(eps=thresh,zk=zk,sources=sources,charges=charges,targets=targets,pgt=1)
 vals.pottarg

src/FMM2D.jl

@@ -9,8 +9,9 @@ All N-body codes return output in an `FMMVals` structure.
 See documentation of N-body codes for details.
 
 N-body interactions with the Helmholtz kernel
-- [`hfmm3d`](@ref): ``O(N)`` fast mutlipole code
+- [`hfmm3d`](@ref): ``O(N)`` fast mutlipole code for the Helmholtz kernel
 - [`h3ddir`](@ref): ``O(N^2)`` direct code
+- [`lfmm3d`](@ref): ``O(N)`` fast mutlipole code for the Laplace kernel
 
 """
 module FMM2D
@@ -20,6 +21,9 @@ using FMM2D_jll
 
 # Exporting interface
 export hfmm2d
+export lfmm2d, rfmm2d, cfmm2d
+# export bhfmm
+# export stfmm2d
 
 # Fortran input/return types
 Fd = Ref{Float64}
@@ -29,6 +33,7 @@ Fc = Ref{ComplexF64}
 # Common input types
 TFN = Union{Array{Float64},Nothing}
 TCN = Union{Array{ComplexF64},Nothing}
+TFCN = Union{Array{Float64},Array{ComplexF64},Nothing}
 
 # Return struct
 mutable struct FMMVals
@@ -48,5 +53,6 @@ end
 include("helper_functions.jl")
 # Include wrappers
 include("helmholtz_wrappers.jl")
+include("laplace_wrappers.jl")
 
 end

src/helmholtz_wrappers.jl

@@ -53,7 +53,7 @@ in the order: ``\\partial_{xx}``, ``\\partial_{yy}``, ``\\partial_{xy}``, ``\\pa
 Non-zero values may indicate insufficient memory available. See the documentation for the FMM2D library.
 If not set (`nothing`), then FMM2D library was never called.
 """
-function hfmm2d(eps::Float64,zk::Union{Float64,ComplexF64},sources::Array{Float64};
+function hfmm2d(;eps::Float64,zk::Union{Float64,ComplexF64},sources::Array{Float64},
         charges::TCN=nothing,dipvecs::TFN=nothing,dipstr::TCN=nothing,targets::TFN=nothing,
         pg::Integer=0,pgt::Integer=0,nd::Integer=1)
 
@@ -110,9 +110,9 @@ function hfmm2d(eps::Float64,zk::Union{Float64,ComplexF64},sources::Array{Float6
 
     if pg == 3
         if nd > 1
-            hess = zeros(ComplexF64,nd,4,n)
+            hess = zeros(ComplexF64,nd,3,n)
         else
-            hess = zeros(ComplexF64,4,n)
+            hess = zeros(ComplexF64,3,n)
         end
     end
 
@@ -134,9 +134,9 @@ function hfmm2d(eps::Float64,zk::Union{Float64,ComplexF64},sources::Array{Float6
 
     if pgt == 3
         if nd > 1
-            hesstarg = zeros(ComplexF64,nd,4,nt)
+            hesstarg = zeros(ComplexF64,nd,3,nt)
         else
-            hesstarg = zeros(ComplexF64,4,nt)
+            hesstarg = zeros(ComplexF64,3,nt)
         end
     end
 
@@ -171,8 +171,7 @@ end
 
 """
 ```julia
-    vals = h2ddir(zk,sources,targets;
-                    charges=nothing,dipstr=nothing,dipvecs=nothing,pgt=0,nd=1,thresh=1e-16)
+    vals = h2ddir(zk,sources,targets, charges=nothing,dipstr=nothing,dipvecs=nothing,pgt=0,nd=1,thresh=1e-16)
 ```
 This function computes the N-body Helmholtz interactions in two dimensions where the
 interaction kernel is given by ``H_0^{(1)}(kr)`` and its gradients. This is the
@@ -215,7 +214,7 @@ in the order: ``\\partial_{xx}``, ``\\partial_{yy}``, ``\\partial_{xy}``, ``\\pa
 * `vals.gradtarg::Array{ComplexF64}` size (nd,2,nt) or (2,nt) gradient at target locations if requested
 * `vals.hesstarg::Array{ComplexF64}` size (nd,4,nt) or (4,nt) Hessian at target locations if requested
 """
-function h2ddir(zk::Union{ComplexF64,Float64},sources::Array{Float64}, targets::Array{Float64};
+function h2ddir(;zk::Union{ComplexF64,Float64},sources::Array{Float64}, targets::Array{Float64},
                 charges::TCN=nothing,dipvecs::TFN=nothing,dipstr::TCN=nothing,
                 pgt::Integer=0,nd::Integer=1,thresh::Float64=1e-15)
 
@@ -271,9 +270,9 @@ function h2ddir(zk::Union{ComplexF64,Float64},sources::Array{Float64}, targets::
 
     if pgt == 3
         if nd > 1
-            hesstarg = zeros(ComplexF64,nd,4,nt)
+            hesstarg = zeros(ComplexF64,nd,3,nt)
         else
-            hesstarg = zeros(ComplexF64,4,nt)
+            hesstarg = zeros(ComplexF64,3,nt)
         end
     end