Inti

Inti.jl is a Julia library for the numerical solution of boundary and volume integral equations. It offers routines for assembling and solving the linear systems that result from applying the Nyström discretization method. Designed for flexibility and efficiency, the package currently supports the following features:

• Specialized integration routines for computing singular and nearly-singular integrals.
• Integrated support for acceleration routines, including the Fast Multipole Method (FMM) and Hierarchical Matrices, by wrapping external libraries.
• Predefined kernels and integral operators for partial differential equations (PDEs) commonly found in mathematical physics (e.g. Laplace, Helmholtz, Stokes).
• Support for complex geometries in 2D and 3D, either through native parametric representations or by importing mesh files from external sources.
• Efficient construction of complex integral operators from simpler ones through lazy composition.

Installing Julia

Download Julia from julialang.org, or use juliaup installer. We recommend using the latest stable version of Julia, although Inti.jl should work with >=v1.9.

Installing Inti.jl

Inti.jl is registered in the Julia General registry and can be installed by launching a Julia REPL and typing the following command:

]add Inti

Alternatively, one can install the latest version of Inti.jl from the main branch using:

using Pkg; Pkg.add(;url = "https://github.com/IntegralEquations/Inti.jl", rev = "main")

Change rev if a different branch or a specific commit hash is desired.

Installing weak dependencies

Inti.jl comes with a set of optional dependencies that can be installed on demand. These provide additional features which can be useful in certain scenarios (e.g. visualization, meshing, acceleration). For convenience, Inti.jl provides the stack_weakdeps_env! function to install all the weak dependencies at once:

using Inti
Inti.stack_weakdeps_env!(; verbose = false, update = true)

Status `~/.julia/scratchspaces/fb74042b-437e-4c5b-88cf-d4e2beb394d5/weakdeps-1.11/Project.toml`
  [2d63477d] FMM2D v0.2.0
  [1e13804c] FMM3D v1.0.1
  [1a804d9e] FMMLIB2D v0.3.2
  [705231aa] Gmsh v0.3.1
  [8646bddf] HMatrices v0.2.10
  [ee78f7c6] Makie v0.21.14
  [eacbb407] Meshes v0.51.22

Note that the first time you run this command, it may take a while to download and compile the dependencies. Subsequent runs will be faster. If preferred, extensions can be manually controlled by Pkg.adding the desired packages from the list above.

Basic usage

Inti.jl can be used to solve a variety of linear partial differential equations by recasting them as integral equations. The general workflow for solving a problem consists of the following steps:

\[ \underbrace{\fbox{Geometry} \rightarrow \fbox{Mesh}}_{\textbf{pre-processing}} \rightarrow \fbox{\color{red}{Solver}} \rightarrow \underbrace{\fbox{Visualization}}_{\textbf{post-processing}}\]

• Geometry: Define the domain of interest using simple shapes (e.g., circles, rectangles) or more complex CAD models.
• Mesh: Create a mesh to approximate the geometry. The mesh is used to define a quadrature and discretize the boundary integral equation.
• Solver: With a mesh and an accompanying quadrature, Inti.jl's routines provide ways to assemble and solve the system of equations arising from the discretization of the integral operators. The core of the library lies in service of this step.
• Visualization: Visualize the solution using a plotting library such as Makie.jl, or export it to a file for further analysis.

As a simple example illustrating the steps above, consider an interior Laplace problem, in two dimensions, with Dirichlet boundary conditions:

\[\begin{aligned} \Delta u &= 0 \quad \text{in } \Omega,\\ u &= g \quad \text{on } \Gamma, \end{aligned}\]

where $\Omega \subset \mathbb{R}^2$ is a sufficiently smooth domain, and $\Gamma = \partial \Omega$ its boundary. A boundary integral reformulation can be achieved by e.g. searching for the solution $u$ in the form of a single-layer potential:

\[u(\boldsymbol{r}) = \int_\Gamma G(\boldsymbol{r},\boldsymbol{y})\sigma(\boldsymbol{y}) \ \mathrm{d}\Gamma(\boldsymbol{y}),\]

where $\sigma : \Gamma \to \mathbb{R}$ is an unknown density function, and $G$ is the fundamental solution of the Laplace equation. This ansatz is, by construction, an exact solution to the PDE on $\Omega$. Imposing the boundary condition on $\Gamma$ leads to the following integral equation:

\[ \int_\Gamma G(\boldsymbol{x},\boldsymbol{y})\sigma(\boldsymbol{y}) \ \mathrm{d}\Gamma(\boldsymbol{y}) = g(\boldsymbol{x}), \quad \forall \boldsymbol{x} \in \Gamma.\]

Expressing the problem above in Inti.jl looks like this:

using Inti, LinearAlgebra, StaticArrays
# create a geometry given by a function f : [0,1] → Γ ⊂ R^2.
geo = Inti.parametric_curve(0, 1) do s
    SVector(0.25, 0.0) + SVector(cos(2π * s) + 0.65 * cos(4π * s[1]) - 0.65, 1.5 * sin(2π * s))
end
Γ = Inti.Domain(geo)
# create a mesh and quadrature
msh = Inti.meshgen(Γ; meshsize = 0.1)
Q = Inti.Quadrature(msh; qorder = 5)
# create the integral operators
op = Inti.Laplace(;dim=2)
S, _ = Inti.single_double_layer(;
    op,
    target = Q,
    source = Q,
    compression = (method = :none,),
    correction = (method = :dim,)
)
# manufacture a harmonic function (exact solution) and take its trace on Γ
uₑ = x -> x[1] + x[2] + x[1]*x[2] + x[1]^2 - x[2]^2  - 2 * log(norm(x .- SVector(-0.5, -1.5)))
g = map(q -> uₑ(q.coords), Q) # value at quad nodes
# solve for σ
σ = S \ g
# use the single-layer potential to evaluate the solution
𝒮, 𝒟 = Inti.single_double_layer_potential(; op, source = Q)
uₕ = x -> 𝒮[σ](x)

#7 (generic function with 1 method)

The function uₕ is now a numerical approximation of the solution to the Laplace equation, and can be evaluated at any point in the domain:

pt = SVector(0.5, 0.1)
println("Exact value at $pt:   ", uₑ(pt))
println("Approx. value at $pt: ", uₕ(pt))

Exact value at [0.5, 0.1]:   -0.3797605448639393
Approx. value at [0.5, 0.1]: -0.3797604679893083

If we care about the solution on the entire domain, we can visualize it using:

using Meshes, GLMakie # trigger the loading of some Inti extensions
xx = yy = range(-2, 2, length = 100)
fig = Figure(; size = (600,300))
inside = x -> Inti.isinside(x, Q)
opts = (xlabel = "x", ylabel = "y", aspect = DataAspect())
ax1 = Axis(fig[1, 1]; title = "Exact solution", opts...)
h1 = heatmap!(ax1, xx,yy,(x, y) -> inside((x,y)) ? uₑ((x,y)) : NaN)
viz!(msh; segmentsize = 3)
cb = Colorbar(fig[1, 3], h1, size = 20, height = 200)
ax2 = Axis(fig[1, 2]; title = "Approx. solution", opts...)
h2 = heatmap!(ax2, xx,yy, (x, y) -> inside((x,y)) ? uₕ((x,y)) : NaN, colorrange = cb.limits[])
viz!(msh; segmentsize = 3)

While the example above is a simple one, Inti.jl can handle significantly more complex problems involving multiple domains, heterogeneous coefficients, vector-valued PDEs, and three-dimensional geometries. The best way to dive deeper into Inti.jl's capabilities is the tutorials section. More advanced usage can be found in the examples section.

Contributing

There are several ways to contribute to Inti.jl:

• Reporting bugs: If you encounter a bug, please open an issue on the GitHub. If possible, please include a minimal working example that reproduces the problem.
• Examples: If you have a cool example that showcases Inti.jl's capabilities, consider submitting a PR to add it to the examples section.
• Contributing code: If you would like to contribute code to Inti.jl, please fork the repository and submit a pull request. Feel free to open a draft PR early in the development process to get feedback on your changes.
• Feature requests: If you have an idea for a new feature or improvement, we would love to hear about it.
• Documentation: If you find any part of the documentation unclear or incomplete, please let us know. Or even better, submit a PR with the improved documentation.

Acknowledgements