Integration of functions over geometries

Specific Domains
, ,
1

Cross-posting from Zulip’s #meshes.jl channel…

Does anyone know of a good reference for integration rules over different geometries?

For example, to integrate a function f over a triangle, we can consider the surface integral:

\int_S fdS = \int_{u=0}^1 \int_{v=0}^{1-u} f(u,v) ||r_u \times r_v|| dudv = Area\ \frac{f_1 + f_2 + f_3}{3}
​
where r = (x(u,v), y(u,v), z(u,v)) is a point in the triangle expressed in terms of parametric coordinates u and v.

More specifically, I am looking for formulas for bilinear integration over quadrangles and trilinear integration over hexahedron in terms of the vertex values f_1, f_2, f_3, ... and their coordinates x_1, x_2, ....

1 Like
2

P.C. Hammer, A.H. Stroud, Numerical integration over simplexes. Math. Tables Aids
Comput. 10 , 137–139 (1956).

M. Abramowitz, I. Stegun, Handbook of Mathematical Functions , 9th edn. (Dover Pub-
lications Inc., New York, NY, 1972).

P.J. Davis, P. Rabinowitz, Methods of Numerical Integration . Computer Science and
Applied Mathematics (Academic Press, New York, NY, 1975).

H. Brass, K. Petras, Quadrature Theory , vol. 178, Mathematical Surveys and Monographs
(American Mathematical Society, Providence, RI, 2011).

R. Cools, P. Rabinowitz, Monomial cubature rules since “Stroud”: a compilation. J. Com-
put. Appl. Math. 48 (3), 309–326 (1993).

R. Cools, Monomial cubature rules since “Stroud”: a compilation. II. J. Comput. Appl.
Math. 112 (1-2), 21–27 (1999).

D.A. Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the tri-
angle. Internat. J. Numer. Methods Engrg. 21 (6), 1129–1148 (1985).

P.C. Hammer, A.H. Stroud, Numerical integration over simplexes. Math. Tables Aids
Comput. 10 , 137–139 (1956).

A.H. Stroud, Approximate Calculation of Multiple Integrals (Prentice-Hall Inc., Engle-
wood Cliffs, NJ, 1971).

M. Abramowitz, I. Stegun, Handbook of Mathematical Functions , 9th edn. (Dover Pub-
lications Inc., New York, NY, 1972).

6 Likes
3

has a medley of (relatively) low-order rules. But not Gauss-Lobatto (useful for nodal integration in spectral methods), except for the lowest order.

4

(a native Julia version of this would be great…)

1 Like
5

GitHub - JuliaApproximation/FastGaussQuadrature.jl: Julia package for Gaussian quadrature?

6

Amazing links. Thank you all. The plan is to implement some of these for Meshes.jl geometries.

3 Likes
7

@Perrin_Meyer is quadpy closed-source? Can’t find their source code to understand what is going on.

8

@Oscar_Smith I think FastGaussQuadrature.jl doesn’t implement integration over geometries. It is mostly hypercubes as far as I know.

9

This is a very good C++ library for integration over convex polytopes that I used some time ago. No Julia interface as far as I know, however…

1 Like
10

You can take a quadrangle or a hexahedron and easily map into into a hypercube by a coordinate transformation, at which point you can apply any quadrature rule for a hypercube (with a corresponding Jacobian factor), whether it is a simple first-order bilinear rule or a tensor product of Gaussian quadrature rules or an adaptive quadrature scheme.

More generally, things like triangular domains can also be converted to rectangles with a change of variables, and general polygons could be decomposed into triangles, etcetera.

2 Likes
11

Weird, because last time I looked (six months ago) all the python code was right on github, and I thought it was all open source… My best guess is this is some sort of code-reorganization?

12

Yes, it is unfortunate that the source code is not open, it would be a great resource to understand the tradeoff of different methods and designs.

13

@juliohm Did you ever come around to making an integration scheme work on Meshes.jl? I can’t seem to find anything…

14

@RomeoV we added some basic functionality in the GeoStats.jl modules on top of Meshes.jl. I think there is a integrate utility function in GeoStatsBase.jl if I remember correctly.

There is also some ongoing work by @mikeingold on line integrals with unitful coordinates over 1D polytopes from Meshes.jl: https://github.com/mikeingold/LineIntegrals.jl

We are in the process of generalizing the coordinate infrastructure in Meshes.jl to accommodate units and non Euclidean coordinates.

1 Like
15

As-is anyone interested in playing around with it can Pkg.add LineIntegals directly from GitHub. The core functionality is to link up Meshes.jl with QuadGK.jl to enable line/path/contour integrals where the geometry is specified using Meshes, e.g.:

using Meshes
using QuadGK
using LineIntegrals

# Construct a path that approximates a unit circle on the xy-plane
#   embedded in a 3D coordinate system
unit_circle = BezierCurve(
    [Point(cos(t), sin(t), 0.0) for t in range(0,2pi,length=361)]
)

# Real function
f(x, y, z) = abs(x + y)
f(p::Point) = f(p.coords...)
quadgk(f, unit_circle)       # -> (5.551055333711397, 1.1102230246251565e-16)

The implementation for 1D polytopes is actually very simple since Meshes already defines a parametric function (::Segment)(t) -> ::Point. The new methods just define a parameterized 1D integral from 0 to 1, and as long as you have a method f(::Point) defined it just works.

I’ve done some testing and everything seems to also compose well with dimensionful types like those from Unitful.jl and DynamicQuantities.jl so you could also have, say, a f that outputs values in units Ohms/meter, integrate f along some Meshes geometry with dimensionful coordinates and the quadgk result will be dimensionful in units Ohms.

I’d considered making this a standalone utility package for General, but am now leaning toward just trying to upstream the useful bits into a Meshes extension. Defining QuadGK.quadgk(f, ::Meshes.*) from a third-party package is type-piracy; I could technically work around this by just defining some LineIntegrals.integral(f, ::Meshes.*) but I don’t know that I’d be adding any real value with a thin rebranded wrapper. Probably preferable to just have a Meshes extension defined for QuadGK that enables the example above without a third package dependency.

1 Like
16

Thanks @mike.ingold for leading this effort. Have you considered

as an alternative to QuadGK.jl? It is native Julia and is super fast from what I remember.

17

I’ve heard of it but haven’t yet had a chance to try it out.

I’ve tended towards QuadGK for 1D problems in the past largely because it’s native Julia that seems to compose well with arbitrary f output types (basically anything that supports (+,-,*,norm) operations) and “just works” for arbitrary atol or rtol specs. The downside I’ve noticed with it, and other adaptive quadrature/cubature implementations, is that if you treat them like a magic bullet and happen to have an f with certain behaviors (like a step change) then solution times will go through the roof.

I see that the FastGaussQuadrature.jl docs promise “16-digit accuracy and in O(n) time”, which seems pretty promising. It looks like all of the built-in quadratures operate on fixed domains, so I’d probably need to implement a change of variables internally to map t instead to [-1,1].

I’ll definitely check it out. Thanks!

1 Like
18

Also, thanks to @RomeoV for bumping this thread I didn’t realize existed. Longer term, I’m also very interested in having a way to compute surface integrals over arbitrary meshes and 2D geometries. I’ll definitely be reading further into the references SGJ provided above.

1 Like
19

I did some more research into FastGaussQuadrature.jl yesterday; the paper they reference was a nice read to help me get up to speed. It looks like the implicit caveat to the “16-digit accuracy and in O(n) time” bumper sticker is that this is only precisely true for functions that can be exactly represented by a polynomial of order 2n-1 over the specified interval. Off the top of my head, I’d guess that for functions with periodic content you’d want to configure n to a bare minimum of the number of expected periods per line/geometric feature, i.e, n>l/\lambda.

One of the features I really like about adaptive quadratures like QuadGK.jl is that a quadgk call can be made without any explicit keyword arguments and the solver will produce a result with some “guaranteed” error tolerance, and also explicitly provide an estimate of the actual resultant error tolerance. Granted, from what I’ve now learned, it sounds like these error estimates are by default predicated on the assumption that the given f is well-approximated over the specified interval by a polynomial of order 15 (order=7).

Within LineIntegrals.jl I originally had methods integral(f, geometry) which were really just wrappers over quadgk that returned only the integral value. That seemed like a weak abstraction since it was really just passing through kwargs and only picking off the result[1], so I re-implemented those as simple QuadGK.quadgk methods.

That being said, I definitely think a case could be made that implementing an interface integral(f, geometry) that internally applies tools from FastGaussQuadrature.jl is a non-trivial value-add. I only had a few minutes last night, but I made a quick-and-dirty integral(f, ::Segment) method that works for scalar functions; I’m going to take another pass to generalize it and make sure I’m minimizing allocations (maybe later implement an in-place integrate! version). The only thing I’m not sure about is how to set a default order n. I’m leaning towards picking a reasonable/conservative default like other packages do (7?) and just making sure that the decision criteria is clear in the docstring.

1 Like
20

That is very exciting @mike.ingold. This approach where you define the high-level interface as integral(f, g) and workout the details with different backends is great. I would suggest that backend options are passed as kwargs so that users can fine-tune the results in specific applications. Something like integral(f, g, order = 7).

If QuadGK.jl has comparable computational performance and better numerical performance, then it makes sense to use it as the default Julia native backend.

Looking forward to seeing this effort evolve into a whole set of integrals over Meshes.jl geometries :100: All our geometries are already parametric in the sense that you can g(u, v) and g(u, v, w) with unitless coordinates u, v, w in [0,1].

1 Like