# Numerical integration over 3D domain of vector-valued function

**URL:** <https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284>\
**Category:** General Usage\
**Tags:** integral\
**Created:** [March 29, 2024, 12:27pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284 "2024-03-29T12:27:23Z")\
**Posts on this page:** 14\
**Page:** 2

<div class="post-metadata">

**Author:** ![mike.ingold](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mike.ingold/32/203749_2.png) [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)\
**Post date:** [April 2, 2024, 4:27pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/21 "2024-04-02T16:27:25Z")

</div>

> [@Domenico\_Lahaye](#):
>
> Sincere thanks! I wonder:
> 
> 1/ Does [Meshes.jl](https://juliahub.com/ui/Packages/General/Meshes) allow to encapsulated meshes generated by [GMSH.jl](https://juliahub.com/ui/Packages/General/GMSH) ? Are there examples of such a use case?
> 
> 2/ Can [Meshes.jl](https://juliahub.com/ui/Packages/General/Meshes) be instructed to recognize regions in the mesh that [GMSH.jl](https://juliahub.com/ui/Packages/General/GMSH) creates?

~~@juliohm is probably better positioned to answer these questions.~~ Edit: he beat me to it lol.

> [@Domenico\_Lahaye](#):
>
> 3/ Can [MeshIntegrals.jl](https://juliahub.com/ui/Packages/General/MeshIntegrals) be extended to treat a singular part of the integral analytically?
> 
> 4/ Can [MeshIntegrals.jl](https://juliahub.com/ui/Packages/General/MeshIntegrals) be extended to treat integrals that depend on a parameter allowing these integrals to be differentiated in a next step?

MeshIntegrals.jl is a fairly new package (started earlier this year). The current interface is simple, i.e. `integral(f, ::Meshes.Geometry, ::IntegrationAlgorithm)`, where you simply pass in some function with a method `f(::Meshes.Point)`. Feel free to submit a PR or an Issue for specific feature requests

Currently:

- If your function has singularities on the integration boundaries then the underlying integration libraries, like QuadGK.jl, should still be able to handle it. If you have singularities within the integration domain then you’d currently need to handle that externally.
- If you have some parameters that are spatially-invariant, then (currently) you’d have to apply them externally, e.g. something like

```julia
f(pt::Point, params) = ...
integral(pt -> f(pt, myparams), geometry, alg)

```

- If your parameters are spatially-variant then you could do the same but with a wrapper function like `myparams(pt::Point)` that returns the appropriate set of parameters based on the current Point.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [April 2, 2024, 4:30pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/22 "2024-04-02T16:30:32Z")

</div>

Many thanks to both @juliohm and @mike.ingold . I will study the matter and get back to you!

---

<div class="post-metadata">

**Author:** ![kylebeggs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kylebeggs/32/43348_2.png) [@kylebeggs](https://discourse.julialang.org/u/kylebeggs)\
**Post date:** [April 4, 2024, 1:19pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/23 "2024-04-04T13:19:25Z")

</div>

If you know where the singularities are, you can simply break the integration into multiple separate integrals where the singularities are the boundaries of each and then sum the results.

---

<div class="post-metadata">

**Author:** ![lxvm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lxvm/32/50010_2.png) [@lxvm](https://discourse.julialang.org/u/lxvm)\
**Post date:** [April 4, 2024, 2:49pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/24 "2024-04-04T14:49:02Z")

</div>

An adaptive integrator given a function with singularities at the integration endpoints can be made much more efficient if the singularity can be subtracted from the function, so the non-singular difference can be integrated numerically, and then its analytic integral added back. The [QuadGK.jl manual](https://juliamath.github.io/QuadGK.jl/stable/quadgk-examples/#Cauchy-principal-values) has some examples of this.

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [April 4, 2024, 8:03pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/25 "2024-04-04T20:03:24Z")

</div>

If I got it right, you are trying to integrate the function  
\int\_V \frac{r}{||r - \tilde{r}||} d\tilde{r}, which is basically r times a convolution of Laplace Green’s Function with a step function. With that insight there are extremely simple things to dom (and quite efficient).

First, compute the convolution separately, in that way you compute 1 difficult integral instead of 3 integrals.

Now, note that \Delta \frac{x^2 + y^2 + z^2 }{6} = 1, so you have that  
\int\_V \frac{1}{||r - \tilde{r}||} d\tilde{r} = \frac{1}{6}\int\_V \frac{\Delta\_{\tilde{r}} \tilde{r}^2}{||r - \tilde{r}||} d\tilde{r}.

By using Green’s third identity, you can reduce the integral to a surface integral problem, i.e.

\frac{1}{6}\int\_V \frac{\Delta\_{\tilde{r}} \tilde{r}^2}{||r - \tilde{r}||} d\tilde{r} = -\frac{r^2}{6} + \int\_{\Gamma} \frac{1}{||r - \tilde{r}||} \frac{\tilde{r}\cdot dS}{3} - \int\_{\Gamma} \frac{\tilde{r}^2}{6} \frac{r - \tilde{r}}{||r - \tilde{r}||^3} \cdot dS

Since this is a cube the problem is reduced to compute 6 surface integral for each r.  
Note that they only behave badly when you get to close to the surface, several things can be done in this case, some are nicer than others, but I would start just using quadgk.  
(Also: if you are interested just in the case for r \in \partial V aka the boundary of the cube, you can do more stuff noticing that \Delta {\bf c} \cdot \tilde{r} = 0 for all \bf c)

(Idea from [[2209.03844] Fast, high-order numerical evaluation of volume potentials via polynomial density interpolation](https://arxiv.org/abs/2209.03844))

Second option, if you don’t want to many digits and you are only interested in values far from the boundary of the cube, you can just use the fact that you have a convolution in free space. Take Fourier transform of both quantities, multiply, and then take the inverse Fourier transform. Some care has to be taken, because the fourier transform of 1/||x|| explodes at 0. Basically follow the steps of this paper [[1604.03155] Fast convolution with free-space Green's functions](https://arxiv.org/abs/1604.03155)

This method, though fast, introduces aliasing near the faces of the cube, so its way less accurate than the previous one.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [April 5, 2024, 7:09pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/26 "2024-04-05T19:09:54Z")

</div>

Correct. Other examples of singularity extraction are in the paper

Duffy, M. G., “Quadrature over a pyramid or cube of integrands with a singularity at a vertex,” SIAM Journal on Numerical Analysis, Vol. 19, 1260–1262, December 1982.

and 500+ papers that cite Duffy-1982.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [April 5, 2024, 7:13pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/27 "2024-04-05T19:13:17Z")

</div>

> [@Veenty](#):
>
> but I would start just using quadgk

I wonder: how do you see quakgk fit to integrate over 2D sufaces? I might miss the obvious here. Thx.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [April 5, 2024, 7:18pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/28 "2024-04-05T19:18:26Z")

</div>

> [@kylebeggs](#):
>
> you can simply break the integration into multiple separate integrals where the singularities are the boundaries of each and then sum the results.

Thx! I agree. An example is

```julia
# define two input integrand
integrand(x,xp) = (x-xp)/abs(x - xp)^1.5

# compute integral by quadrature over second input
u(x) = quadgk(xp -> integrand(x,xp), 0, 1)

```

u(1) evaluates fine. u(0.9) does not,

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [April 5, 2024, 7:32pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/29 "2024-04-05T19:32:30Z")

</div>

The thing is that as long as r \not \in \partial V, the integrand is a differentiable function, so standard refinement eventually works fine.

But to make it even better, if F\_i are the faces of the cube, for each i, you look for the r^\*\_i \in F\_i that is closest to r, this is where the almost singularity will be. Then divide the face in 4 sections, where r^\* is the point of intersection of this 4 sections. Now you have 4 integrals that they have almost-singularities in one corner.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [April 5, 2024, 8:10pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/30 "2024-04-05T20:10:10Z")

</div>

Thx! Are you suggesting that replacing the volume integral over a cell by a set of surface integrals over the faces allows a better grip on the singularity?

---

<div class="post-metadata">

**Author:** ![Veenty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/veenty/32/50940_2.png) [@Veenty](https://discourse.julialang.org/u/Veenty)\
**Post date:** [April 5, 2024, 8:40pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/31 "2024-04-05T20:40:10Z")

</div>

Yes, see that with the Green’s function trick that I showed, you only have a singularity if the point that you are evaluating is close to the surface of the cube, otherwise the integrand is smooth.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [April 6, 2024, 5:40am UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/32 "2024-04-06T05:40:59Z")

</div>

Clear, thx. If this is the situation, you would like the mesh operations in Julia to support these kind of manipulations (retrieve faces belonging to a cell, subdivide a face, compute distance to nearest vertex etc). Any suggestions on what is already available and how to deploy it?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [April 6, 2024, 11:41am UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/33 "2024-04-06T11:41:02Z")

</div>

All these operations are already supported in Meshes.jl, please check the docs and ask in our Zulip channel if you have specific questions.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [April 6, 2024, 6:52pm UTC](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/34 "2024-04-06T18:52:50Z")

</div>

Thx. Will make homework and get back to you.

[Previous page](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284.md?page=1)