### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [May 31, 2022, 2:04pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/1 "2022-05-31T14:04:05Z")
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, ....
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [May 31, 2022, 2:33pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/2 "2022-05-31T14:33:34Z")
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).
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [May 31, 2022, 2:47pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/3 "2022-05-31T14:47:46Z")
>
has a medley of (relatively) low-order rules. But not Gauss-Lobatto (useful for nodal integration in spectral methods), except for the lowest order.
---
### Author:  [@Perrin\_Meyer](https://discourse.julialang.org/u/Perrin_Meyer)
#### Post date: [May 31, 2022, 4:41pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/4 "2022-05-31T16:41:11Z")
> **[GitHub - sigma-py/quadpy: :triangular\_ruler: Numerical integration...](https://github.com/sigma-py/quadpy)**
>
> :triangular\_ruler: Numerical integration (quadrature, cubature) in Python - GitHub - sigma-py/quadpy: :triangular\_ruler: Numerical integration (quadrature, cubature) in Python
(a native Julia version of this would be great…)
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [May 31, 2022, 4:45pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/5 "2022-05-31T16:45:28Z")
[GitHub - JuliaApproximation/FastGaussQuadrature.jl: Julia package for Gaussian quadrature](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)?
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [May 31, 2022, 4:51pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/6 "2022-05-31T16:51:27Z")
Amazing links. Thank you all. The plan is to implement some of these for Meshes.jl geometries.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [May 31, 2022, 5:48pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/7 "2022-05-31T17:48:11Z")
@Perrin_Meyer is quadpy closed-source? Can’t find their source code to understand what is going on.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [May 31, 2022, 5:48pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/8 "2022-05-31T17:48:46Z")
@Oscar_Smith I think FastGaussQuadrature.jl doesn’t implement integration over geometries. It is mostly hypercubes as far as I know.
---
### Author:  [@Arrigo\_Benedetti](https://discourse.julialang.org/u/Arrigo_Benedetti)
#### Post date: [May 31, 2022, 5:49pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/9 "2022-05-31T17:49:04Z")
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…
> **[GitHub - GeomScale/volesti: Practical volume computation and sampling in high...](https://github.com/GeomScale/volesti)**
>
> Practical volume computation and sampling in high dimensions - GitHub - GeomScale/volesti: Practical volume computation and sampling in high dimensions
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 31, 2022, 7:05pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/10 "2022-05-31T19:05:09Z")
> [@juliohm](#):
>
> 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, \ldots and their coordinates x\_1, x\_2, \ldots .
You can take a quadrangle or a hexahedron and easily map into into a hypercube by a [coordinate transformation](https://www.geometrictools.com/Documentation/PerspectiveMappings.pdf), 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](https://github.com/JuliaMath/HCubature.jl) scheme.
More generally, things like triangular domains can also be [converted to rectangles](https://discourse.julialang.org/t/2d-integration-over-non-rectangular-domain-using-cubature/2991/7) with a change of variables, and general polygons could be [decomposed into triangles](https://discourse.julialang.org/t/integral-of-a-function-inside-non-rectangular-domain/72132/9), etcetera.
---
### Author:  [@Perrin\_Meyer](https://discourse.julialang.org/u/Perrin_Meyer)
#### Post date: [May 31, 2022, 9:50pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/11 "2022-05-31T21:50:31Z")
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?
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [June 1, 2022, 3:38pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/12 "2022-06-01T15:38:32Z")
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.
---
### Author:  [@RomeoV](https://discourse.julialang.org/u/RomeoV)
#### Post date: [January 29, 2024, 8:46am UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/13 "2024-01-29T08:46:30Z")
@juliohm Did you ever come around to making an integration scheme work on Meshes.jl? I can’t seem to find anything…
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [January 29, 2024, 9:09am UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/14 "2024-01-29T09:09:09Z")
@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](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.
---
### Author:  [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)
#### Post date: [January 29, 2024, 2:10pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/15 "2024-01-29T14:10:14Z")
> [@juliohm](#):
>
> There is also some ongoing work by @mikeingold on line integrals with unitful coordinates over 1D polytopes from [Meshes.jl](https://juliahub.com/ui/Packages/Meshes): [https://github.com/mikeingold/LineIntegrals.jl](https://github.com/mikeingold/LineIntegrals.jl)
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.:
```julia
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.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [January 29, 2024, 2:31pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/16 "2024-01-29T14:31:17Z")
Thanks @mike.ingold for leading this effort. Have you considered
> **[Home · FastGaussQuadrature.jl](https://juliaapproximation.github.io/FastGaussQuadrature.jl/stable/)**
>
> Documentation for FastGaussQuadrature.jl.
as an alternative to QuadGK.jl? It is native Julia and is super fast from what I remember.
---
### Author:  [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)
#### Post date: [January 29, 2024, 3:05pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/17 "2024-01-29T15:05:57Z")
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!
---
### Author:  [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)
#### Post date: [January 29, 2024, 3:20pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/18 "2024-01-29T15:20:12Z")
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.
---
### Author:  [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)
#### Post date: [January 30, 2024, 2:43pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/19 "2024-01-30T14:43:40Z")
I did some more research into FastGaussQuadrature.jl yesterday; the [paper they reference](http://pi.math.cornell.edu/~ajt/papers/QuadratureEssay.pdf) 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](https://github.com/mikeingold/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.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [January 30, 2024, 3:01pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/20 "2024-01-30T15:01:05Z")
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 💯 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]`.
[Next page](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994.md?page=2)