Skip to content

Commit ab731fd

Browse files
maltezfariamaltezfaria
committed
add back singular quadrature rules
1 parent fe1fe3a commit ab731fd

5 files changed

Lines changed: 406 additions & 0 deletions

File tree

‎src/Integration/Integration.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,3 +7,5 @@ singular integration routines useful for (weakly) singular integrands.
77

88
include("quadrulestables.jl")
99
include("quadrule.jl")
10+
include("singularityhandler.jl")
11+
include("singularquadrule.jl")

‎src/Integration/singularityhandler.jl‎

Lines changed: 197 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,197 @@
1+
"""
2+
abstract type AbstractSingularityHandler{R}
3+
4+
Used for handling localized integrand singularities in `R`.
5+
"""
6+
abstract type AbstractSingularityHandler{T} end
7+
8+
"""
9+
struct IMT{A,P} <: AbstractSingularityHandler{ReferenceLine}
10+
11+
One-dimensional change of variables mapping `[0,1] -> [0,1]` with the property that
12+
all derivatives vanish at the point `x=0`.
13+
14+
See [Davis and Rabinowitz](https://www.elsevier.com/books/methods-of-numerical-integration/davis/978-0-12-206360-2).
15+
"""
16+
struct IMT{A,P} <: AbstractSingularityHandler{ReferenceLine}
17+
end
18+
IMT(;a=1,p=1) = IMT{a,p}()
19+
20+
domain(::IMT) = ReferenceLine()
21+
22+
"""
23+
image(f)
24+
25+
The a function-like object `f: Ω → R`, return `R`.
26+
"""
27+
image(::IMT) = ReferenceLine()
28+
29+
function (f::IMT{A,P})(x) where {A,P}
30+
exp(A * (1 - 1 / x[1]^P))
31+
end
32+
33+
derivative(f::IMT{A,P},x) where {A,P} = f(x) * A * P * 1 / x[1]^(P + 1)
34+
35+
jacobian(f::IMT,x) = derivative(f, x) |> SMatrix{1,1}
36+
37+
"""
38+
struct Kress{P} <: AbstractSingularityHandler{ReferenceLine}
39+
40+
Change of variables mapping `[0,1]` to `[0,1]` with the property that the first
41+
`P-1` derivatives of the transformation vanish at `x=0`.
42+
"""
43+
struct Kress{P} <: AbstractSingularityHandler{ReferenceLine}
44+
end
45+
Kress(;order=5) = Kress{order}()
46+
47+
domain(k::Kress) = ReferenceLine()
48+
image(k::Kress) = ReferenceLine()
49+
50+
# NOTE: fastmath is needed here to allow for various algebraic simplifications
51+
# which are not exact in floating arithmetic. Maybe reorder the operations *by
52+
# hand* to avoid having to use fastmath? In any case, benchmark first.
53+
function (f::Kress{P})(y) where {P}
54+
x = y[1]
55+
v = (x) -> (1 / P - 1 / 2) * ((1 - x))^3 + 1 / P * ((x - 1)) + 1 / 2
56+
return 2v(x)^P / (v(x)^P + v(2 - x)^P)
57+
end
58+
59+
function derivative(f::Kress{P}, y) where {P}
60+
x = y[1]
61+
v = (x) -> (1 / P - 1 / 2) * ((1 - x))^3 + 1 / P * ((x - 1)) + 1 / 2
62+
vp = (x) -> -3 * (1 / P - 1 / 2) * ((1 - x))^2 + 1 / P
63+
return 2 * (P * v(x)^(P - 1) * vp(x) * (v(x)^P + v(2 - x)^P) - (P * v(x)^(P - 1) * vp(x) - P * v(2 - x)^(P - 1) * vp(2 - x) ) * v(x)^P ) /
64+
(v(x)^P + v(2 - x)^P)^2
65+
end
66+
67+
jacobian(f::Kress,x) = derivative(f, x) |> SMatrix{1,1}
68+
69+
struct KressR{P} <: AbstractSingularityHandler{ReferenceLine}
70+
end
71+
KressR(;order=5) = KressR{order}()
72+
73+
domain(k::KressR) = ReferenceLine()
74+
image(k::KressR) = ReferenceLine()
75+
76+
# NOTE: fastmath is needed here to allow for various algebraic simplifications
77+
# which are not exact in floating arithmetic. Maybe reorder the operations *by
78+
# hand* to avoid having to use fastmath? In any case, benchmark first.
79+
function (f::KressR{P})(y) where {P}
80+
x = 1-y[1]
81+
v = (x) -> (1 / P - 1 / 2) * ((1 - x))^3 + 1 / P * ((x - 1)) + 1 / 2
82+
return 1 - 2v(x)^P / (v(x)^P + v(2 - x)^P)
83+
end
84+
85+
function derivative(f::KressR{P}, y) where {P}
86+
x = 1-y[1]
87+
v = (x) -> (1 / P - 1 / 2) * ((1 - x))^3 + 1 / P * ((x - 1)) + 1 / 2
88+
vp = (x) -> -3 * (1 / P - 1 / 2) * ((1 - x))^2 + 1 / P
89+
return 2 * (P * v(x)^(P - 1) * vp(x) * (v(x)^P + v(2 - x)^P) - (P * v(x)^(P - 1) * vp(x) - P * v(2 - x)^(P - 1) * vp(2 - x) ) * v(x)^P ) /
90+
(v(x)^P + v(2 - x)^P)^2
91+
end
92+
93+
jacobian(f::KressR,x) = derivative(f, x) |> SMatrix{1,1}
94+
95+
"""
96+
struct KressP{P} <: AbstractSingularityHandler{ReferenceLine}
97+
98+
Like [`Kress`](@ref), this change of variables maps the interval `[0,1]` onto
99+
itself, but the first `P` derivatives of the transformation vanish at **both**
100+
endpoints.
101+
102+
This change of variables can be used to *periodize* integrals in the following
103+
sense. Suppose we wish to compute the integral of `f(x)` from `0` to `1` where
104+
`f is not a `1-periodic function. If `ϕ` is an object of type `KressP`, then
105+
using it as a change of variables in the integration yields a similar integral
106+
from `0` to `1` (the interval `0≤0≤1` is mappend onto itself), but with
107+
integrand given by `g(x) = f(ϕ(x))ϕ'(x)`. Since `ϕ'` vanishes (together with `P`
108+
of its derivatives), the function `g(x)` is now periodic (up to derivatives of
109+
order up to `P`) at the endpoints. Thus quadrature rules designed for periodic
110+
functions like the [`TrapezoidalOpen`](@ref) can be used to obtain high order
111+
convergence of `g`, which in turn yields a modified quadrature rule when viewed
112+
as a quadrature rule for `f`.
113+
"""
114+
struct KressP{P} <: AbstractSingularityHandler{ReferenceLine}
115+
end
116+
KressP(;order=5) = KressP{order}()
117+
118+
domain(k::KressP) = ReferenceLine()
119+
image(k::KressP) = ReferenceLine()
120+
121+
@fastmath function (f::KressP{P})(y) where {P}
122+
x = y[1]
123+
v = (x) -> (1 / P - 1 / 2) * ((1 - 2x))^3 + 1 / P * ((2x - 1)) + 1 / 2
124+
return v(x)^P / (v(x)^P + v(1 - x)^P)
125+
end
126+
127+
@fastmath function derivative(f::KressP{P}, y) where {P}
128+
x = y[1]
129+
v = (x) -> (1 / P - 1 / 2) * ((1 - 2x))^3 + 1 / P * ((2x - 1)) + 1 / 2
130+
vp = (x) -> -6 * (1 / P - 1 / 2) * ((1 - 2x))^2 + 2 / P
131+
return (P * v(x)^(P - 1) * vp(x) * (v(x)^P + v(1 - x)^P) - (P * v(x)^(P - 1) * vp(x) - P * v(1 - x)^(P - 1) * vp(1 - x) ) * v(x)^P ) /
132+
(v(x)^P + v(1 - x)^P)^2
133+
end
134+
135+
jacobian(f::KressP,x) = derivative(f, x) |> SMatrix{1,1}
136+
137+
"""
138+
struct Duffy <: AbstractSingularityHandler{RefereceTriangle}
139+
140+
Change of variables mapping the `ReferenceSquare` to the `RefereceTriangle` with
141+
the property that the jacobian vanishes at the `(1,0)` vertex of the triangle.
142+
143+
Useful for integrating functions with a singularity on the `(1,0)` edge of the
144+
reference triangle.
145+
"""
146+
struct Duffy <: AbstractSingularityHandler{ReferenceTriangle} end
147+
148+
domain(::Duffy) = ReferenceSquare()
149+
image(::Duffy) = ReferenceTriangle()
150+
151+
function (::Duffy)(u)
152+
SVector(u[1], (1 - u[1]) * u[2])
153+
end
154+
155+
function jacobian(::Duffy, u)
156+
SMatrix{2,2,Float64}(1, 0, -u[2][1], (1 - u[1][1]))
157+
end
158+
159+
# TODO: generalize to `N` dimensions
160+
"""
161+
struct TensorProductSingularityHandler{S} <: AbstractSingularityHandler{ReferenceSquare}
162+
163+
A tensor product of two one-dimensional `AbstractSingularityHandler`s for performing
164+
integration over the `ReferenceSquare`.
165+
"""
166+
struct TensorProductSingularityHandler{S} <: AbstractSingularityHandler{ReferenceSquare}
167+
shandler::S
168+
end
169+
170+
domain(::TensorProductSingularityHandler) = ReferenceSquare()
171+
image(::TensorProductSingularityHandler) = ReferenceSquare()
172+
173+
function TensorProductSingularityHandler(q...)
174+
TensorProductSingularityHandler(q)
175+
end
176+
177+
TensorProductSingularityHandler{Tuple{P,Q}}() where {P,Q} = TensorProductSingularityHandler(P(), Q())
178+
179+
function (f::TensorProductSingularityHandler)(x)
180+
shandler = f.shandler
181+
@assert length(shandler) == length(x)
182+
N = length(shandler)
183+
svector(i -> shandler[i](x[i]), N)
184+
end
185+
186+
function jacobian(f::TensorProductSingularityHandler, x)
187+
shandler = f.shandler
188+
@assert length(shandler) == length(x)
189+
N = length(shandler)
190+
if N == 2
191+
jx = derivative(shandler[1], x[1])
192+
jy = derivative(shandler[2], x[2])
193+
return SMatrix{2,2}(jx, 0, 0, jy)
194+
else
195+
notimplemented()
196+
end
197+
end

‎src/Integration/singularquadrule.jl‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
"""
2+
SingularQuadratureRule{D,Q,S} <: AbstractQuadratureRule{D}
3+
4+
A quadrature rule over `D` intended to integrate functions which are singular at
5+
a known point `s ∈ D`.
6+
7+
A singular quadrature is rule is composed of a *regular* quadrature rule (e.g.
8+
`GaussLegendre`) and a [`AbstractSingularityHandler`](@ref) to transform the regular
9+
quadrature. The regular quadrature rule generates nodes and weights on the
10+
`domain(sing_handler)`, and those are mapped into an appropriate quadrature over
11+
`D = range(sing_handler)` using the singularity handler.
12+
"""
13+
struct SingularQuadratureRule{D,Q,S} <: AbstractQuadratureRule{D}
14+
qrule::Q
15+
singularity_handler::S
16+
function SingularQuadratureRule(q::Q,s::S) where {Q,S}
17+
@assert domain(s) == domain(q) "domain of quadrature must coincide with domain of singularity handler"
18+
D = image(s) |> typeof
19+
new{D,Q,S}(q,s)
20+
end
21+
end
22+
23+
# getters
24+
qrule(q::SingularQuadratureRule) = q.qrule
25+
singularity_handler(q::SingularQuadratureRule) = q.singularity_handler
26+
27+
domain(qs::SingularQuadratureRule) = domain(singularity_handler(qs))
28+
image(qs::SingularQuadratureRule) = image(singularity_handler(qs))
29+
30+
@generated function (qs::SingularQuadratureRule{D,Q,S})() where {D,Q,S}
31+
qstd = Q()
32+
shand = S()
33+
x̂,ŵ = qstd() # reference quadrature
34+
# use the change of variables cov to modify reference quadrature
35+
x = map(x->shand(x),x̂)
36+
w = map(x̂,ŵ) do x̂,ŵ
37+
jac = jacobian(shand,x̂)
38+
μ = abs(det(jac))
39+
μ*prod(ŵ)
40+
end
41+
return :($x,$w)
42+
end
43+
44+
function (qs::SingularQuadratureRule{ReferenceLine})(s)
45+
@assert s ∈ image(qs)
46+
x̂,ŵ = qs()
47+
# left domain
48+
x1 = @. s[1] * (1-x̂)
49+
w1 = @. ŵ*s[1]
50+
# right domain
51+
x2 = @. s[1] + x̂*(1-s[1])
52+
w2 = @. ŵ*(1-s[1])
53+
# combine the nodes and weights
54+
x = vcat(x1,x2)
55+
w = vcat(w1,w2)
56+
return x,w
57+
end
58+
59+
"""
60+
singular_quadrature(k,q::SingularQuadratureRule,s)
61+
62+
Return nodes and weights to integrate a function over `domain(q)` with a
63+
factored weight `k`.
64+
"""
65+
function singular_quadrature(k,q::SingularQuadratureRule,s)
66+
x,w = q(s)
67+
T = Base.promote_op(k,eltype(x))
68+
assert_concrete_type(T)
69+
w = map(zip(x,w)) do (x,w)
70+
k(x)*w
71+
end
72+
return x,w
73+
end
74+
75+
"""
76+
singular_weights(k,xi,q::SingularQuadratureRule,s)
77+
"""
78+
function singular_weights(k,xi,q::SingularQuadratureRule,s)
79+
x,w = singular_quadrature(k,q,s)
80+
x = map(x->x[1],x)
81+
ws = map(lag_basis) do li
82+
integrate(x,w) do x
83+
f(x)*li(x)
84+
end
85+
end
86+
# wlag = barycentric_lagrange_weights(xi)
87+
# L = barycentric_lagrange_matrix(xi,x,wlag)
88+
# return transpose(L)*w
89+
return ws
90+
end
91+
92+
function singular_weights(k,qreg::AbstractQuadratureRule,qsin::SingularQuadratureRule,s)
93+
xq,wq = qsin(s)
94+
xi = qnodes(qreg)
95+
lag_basis = lagrange_basis(xi)
96+
map(lag_basis) do l
97+
integrate(xq,wq) do x
98+
k(x)*l(x)
99+
end
100+
end
101+
# return transpose(L)*w
102+
end

‎test/Integration/runtests.jl‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,3 +4,7 @@ using SafeTestsets
44
@safetestset "Quadrature rules" begin
55
include("quadrule_test.jl")
66
end
7+
8+
@safetestset "Singular quadrature rules" begin
9+
include("singularquadrule_test.jl")
10+
end

0 commit comments

Comments
 (0)