Added DuffyQuadrature

Author: fverdugo

src/Integration/DuffyQuadratures.jl

@@ -0,0 +1,132 @@
+module DuffyQuadratures
+
+using Gridap
+using FastGaussQuadrature: gaussjacobi
+using FastGaussQuadrature: gausslegendre
+using StaticArrays
+
+export DuffyQuadrature
+import Gridap: coordinates
+import Gridap: weights
+
+struct DuffyQuadrature{D} <: Quadrature{D}
+  x::Vector{Point{D,Float64}}
+  w::Vector{Float64}
+end
+
+function DuffyQuadrature{D}(order::Integer) where D
+  x,w = _duffy_quad_data(order,D)
+  DuffyQuadrature{D}(x,w)
+end
+
+coordinates(q::DuffyQuadrature) = q.x
+
+weights(q::DuffyQuadrature) = q.w
+
+# Helpers
+
+function _duffy_quad_data(order::Integer,D::Int)
+
+  beta = 0
+  dim_to_quad_1d = [
+    _gauss_jacobi_in_0_to_1(order,(D-1)-(d-1),beta) for d in 1:(D-1) ]
+
+  quad_1d = _gauss_legendre_in_0_to_1(order)
+  push!(dim_to_quad_1d,quad_1d)
+
+  x_pos = 1
+  w_pos = 2
+  dim_to_xs_1d = [quad_1d[x_pos] for quad_1d in dim_to_quad_1d]
+  dim_to_ws_1d = [quad_1d[w_pos] for quad_1d in dim_to_quad_1d]
+
+  a = 0.5
+  for d in (D-1):-1:1
+    ws_1d = dim_to_ws_1d[d]
+    ws_1d[:] *= a
+    a *= 0.5
+  end
+
+  x,w = _tensor_product(dim_to_xs_1d,dim_to_ws_1d)
+
+  (_duffy_map.(x),w)
+
+end
+
+"""
+Duffy map from the n-cube in [0,1]^d to the n-simplex in [0,1]^d
+"""
+function _duffy_map(q::Point{D,T}) where {D,T}
+  m = zero(MVector{D,T})
+  m[1] = q[1]
+  a = one(T)
+  for i in 2:D
+    a *= (1-q[i-1])
+    m[i] = a*q[i]
+  end
+  Point{D,T}(m)
+end
+
+_duffy_map(q::Point{1,T}) where T = q
+
+function _gauss_jacobi_in_0_to_1(order,alpha,beta)
+  n = _npoints_from_order(order)
+  x,w = gaussjacobi(n,alpha,beta)
+  _map_to(0,1,x,w)
+end
+
+function _gauss_legendre_in_0_to_1(order)
+  n = _npoints_from_order(order)
+  x,w = gausslegendre(n)
+  _map_to(0,1,x,w)
+end
+
+"""
+Transforms a 1-D quadrature from `[-1,1]` to `[a,b]`, with `a<b`.
+"""
+function _map_to(a,b,points,weights)
+  points_ab = 0.5*(b-a)*points .+ 0.5*(a+b)
+  weights_ab = 0.5*(b-a)*weights
+  (points_ab, weights_ab)
+end
+
+function _npoints_from_order(order)
+  ceil(Int, (order + 1.0) / 2.0 )
+end
+
+function _tensor_product(
+  dim_to_xs_1d::Vector{Vector{T}},
+  dim_to_ws_1d::Vector{Vector{W}}) where {T,W}
+
+  D = length(dim_to_ws_1d)
+  @assert D == length(dim_to_xs_1d)
+  dim_to_n = [length(ws_1d) for ws_1d in dim_to_ws_1d]
+  n = prod(dim_to_n)
+  xs = zeros(Point{D,T},n)
+  ws = zeros(W,n)
+  cis = CartesianIndices(tuple(dim_to_n...))
+  m = zero(MVector{D,T})
+  _tensor_product!(xs,ws,dim_to_xs_1d,dim_to_ws_1d,cis,m)
+  (xs,ws)
+end
+
+function _tensor_product!(
+  xs,ws,dim_to_xs_1d,dim_to_ws_1d,cis::CartesianIndices{D},m) where D
+  k = 1
+  for ci in cis
+    w = 1.0
+    for d in 1:D
+      xs_1d = dim_to_xs_1d[d]
+      ws_1d = dim_to_ws_1d[d]
+      i = ci[d]
+      xi = xs_1d[i]
+      wi = ws_1d[i]
+      w *= wi
+      m[d] = xi
+    end
+    xs[k] = m
+    ws[k] = w
+    k += 1
+  end
+end
+
+end # module

src/Integration/files.jl

@@ -2,6 +2,9 @@
 include("Quadratures.jl")
 @reexport using Gridap.Quadratures
 
+include("DuffyQuadratures.jl")
+@reexport using Gridap.DuffyQuadratures
+
 include("CellQuadratures.jl")
 @reexport using Gridap.CellQuadratures
 

test/IntegrationTests/DuffyQuadraturesTests.jl

@@ -0,0 +1,28 @@
+module DuffyQuadraturesTests
+
+using Test
+using Gridap
+
+D = 3
+order=2
+quad = DuffyQuadrature{D}(order)
+@test isa(coordinates(quad),Vector{Point{D,Float64}})
+@test isa(weights(quad),Vector{Float64})
+@test sum(weights(quad)) ≈ 1/6
+
+D = 2
+order=2
+quad = DuffyQuadrature{D}(order)
+x = coordinates(quad)
+r = Point{D,Float64}[
+  (0.155051, 0.178559), (0.644949, 0.0750311),
+  (0.155051, 0.66639), (0.644949, 0.28002)]
+@test isapprox(reinterpret(Float64,x),reinterpret(Float64,r),rtol=1e-3)
+@test sum(weights(quad)) ≈ 1/2
+
+D = 1
+order=2
+quad = DuffyQuadrature{D}(order)
+@test sum(weights(quad)) ≈ 1
+
+end # module

test/IntegrationTests/runtests.jl

@@ -4,6 +4,8 @@ using Test
 
 @testset "Quadratures" begin include("QuadraturesTests.jl") end
 
+@testset "DuffyQuadratures" begin include("DuffyQuadraturesTests.jl") end
+
 @testset "CellQuadratures" begin include("CellQuadraturesTests.jl") end
 
 @testset "Triangulations" begin include("TriangulationsTests.jl") end