Numerical integration over 3D domain of vector-valued function

General Usage
21

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

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
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.
1 Like
22

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

1 Like
23

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.

1 Like
24

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 has some examples of this.

2 Likes
25

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)

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

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

2 Likes
26

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.

27

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

28

Thx! I agree. An example is

# 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,

29

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.

30

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?

31

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.

32

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?

33

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

1 Like
34

Thx. Will make homework and get back to you.

1 Like