Meshes.jl/src/integration.jl at d2e040afdfe42a3f629828e9ef93a6392bd4e3d5 · JuliaGeometry/Meshes.jl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    integral(fun, geom; order=1)

Calculate the integral over the `geom`etry of the `fun`ction that maps
[`Point`](@ref)s to values in a linear space.

The quadrature `order` can be specified to increase the accuracy.

See also [`localintegral`](@ref).
"""
integral(fun, geom; order=1) = localintegral(fun ∘ geom, geom; order)

"""
    localintegral(fun, geom; order=1)

Calculate the integral over the `geom`etry of the `fun`ction that maps
parametric coordinates `uvw` to values in a linear space.

The quadrature `order` can be specified to increase the accuracy.

See also [`integral`](@ref).
"""
function localintegral(fun, geom; order=1)
  # parametric dimension
  dim = paramdim(geom)

  # Gauss-Legendre quadrature points and weights
  ts, ws = gausslegendre(order)
  tgrid = Iterators.product(ntuple(_ -> ts, dim)...)
  wgrid = Iterators.product(ntuple(_ -> ws, dim)...)

  # compute integral with change of variable and differential element
  Σwᵢfᵢ = sum(zip(tgrid, wgrid)) do (t, w)
    # change of variable [-1, 1] → [0, 1]
    uvw = ntuple(i -> (t[i] + 1) / 2, dim)
    prod(w) * fun(uvw...) * differential(geom, uvw)
  end

  # adjust for change of variable
  Σwᵢfᵢ / 2^dim
end