Merge pull request #61 from gridap/duffy_based_integration_for_simplices

Duffy based integration for simplices

Author: fverdugo

Manifest.toml

@@ -3,6 +3,12 @@
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
+[[BinDeps]]
+deps = ["Compat", "Libdl", "SHA", "URIParser"]
+git-tree-sha1 = "12093ca6cdd0ee547c39b1870e0c9c3f154d9ca9"
+uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee"
+version = "0.8.10"
+
 [[BinaryProvider]]
 deps = ["Libdl", "Logging", "SHA"]
 git-tree-sha1 = "c7361ce8a2129f20b0e05a89f7070820cfed6648"
@@ -21,6 +27,12 @@ git-tree-sha1 = "50b3ae4d643dc27eaff69fb6be06ee094d5500c9"
 uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 version = "0.7.0"
 
+[[Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "84aa74986c5b9b898b0d1acaf3258741ee64754f"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "2.1.0"
+
 [[DataStructures]]
 deps = ["InteractiveUtils", "OrderedCollections", "Random", "Serialization", "Test"]
 git-tree-sha1 = "ca971f03e146cf144a9e2f2ce59674f5bf0e8038"
@@ -31,10 +43,20 @@ version = "0.15.0"
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
 [[Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
+[[FastGaussQuadrature]]
+deps = ["Compat", "LinearAlgebra", "SpecialFunctions"]
+git-tree-sha1 = "769ac4057ed875f433372e9a571a2cb86347d1be"
+uuid = "442a2c76-b920-505d-bb47-c5924d526838"
+version = "0.3.3"
+
 [[InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
@@ -62,6 +84,9 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
 [[OrderedCollections]]
 deps = ["Random", "Serialization", "Test"]
 git-tree-sha1 = "c4c13474d23c60d20a67b217f1d7f22a40edf8f1"
@@ -108,13 +133,23 @@ uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
+[[SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
 [[Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
 [[SparseArrays]]
 deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[[SpecialFunctions]]
+deps = ["BinDeps", "BinaryProvider", "Libdl", "Test"]
+git-tree-sha1 = "0b45dc2e45ed77f445617b99ff2adf0f5b0f23ea"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "0.7.2"
+
 [[StaticArrays]]
 deps = ["InteractiveUtils", "LinearAlgebra", "Random", "Statistics", "Test"]
 git-tree-sha1 = "3841b39ed5f047db1162627bf5f80a9cd3e39ae2"
@@ -147,6 +182,12 @@ git-tree-sha1 = "a25d8e5a28c3b1b06d3859f30757d43106791919"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
 version = "0.9.4"
 
+[[URIParser]]
+deps = ["Test", "Unicode"]
+git-tree-sha1 = "6ddf8244220dfda2f17539fa8c9de20d6c575b69"
+uuid = "30578b45-9adc-5946-b283-645ec420af67"
+version = "0.4.0"
+
 [[UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

Project.toml

@@ -5,6 +5,7 @@ version = "0.2.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"

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/QuadratureFactories.jl

@@ -0,0 +1,29 @@
+module QuadratureFactories
+
+using Gridap
+using Gridap.Helpers
+
+import Gridap: Quadrature
+
+"""
+Factory function to create Quadrature objects in a convenient way
+"""
+function Quadrature(p::Polytope,order::Int)
+  Quadrature(p.extrusion.array.data,order=order)
+end
+
+function Quadrature(extrusion::NTuple{D,Int};order::Int) where D
+
+  if all(extrusion .== HEX_AXIS)
+    return TensorProductQuadrature(orders=tuple(fill(order,D)...))
+
+  elseif all(extrusion .== TET_AXIS)
+    return DuffyQuadrature{D}(order)
+
+  else
+    @notimplemented
+  end
+
+end
+
+end # module

src/Integration/Quadratures.jl

@@ -25,15 +25,6 @@ coordinates(::Quadrature{D} where D)::Array{Point{D,Float64},1} = @abstractmetho
 
 weights(::Quadrature)::Array{Float64,1} = @abstractmethod
 
-# Factories
-
-"""
-Factory function to create Quadrature objects in a convenient way
-"""
-function Quadrature(extrusion::NTuple{D,Int};order::Int) where D
-  _quadrature(extrusion,order)
-end
-
 # Concrete structs
 
 """
@@ -62,11 +53,6 @@ weights(self::TensorProductQuadrature) = self.weights
 
 # Helpers
 
-function _quadrature(extrusion::NTuple{D,Int},order) where D
-  @notimplementedif any(extrusion == TET_AXIS)
-  TensorProductQuadrature(orders=tuple(fill(order,D)...))
-end
-
 function _tensor_product(quads,npoints,::Val{D}) where D
   @assert length(quads) == D
   @assert length(npoints) == D

src/Integration/files.jl

@@ -2,6 +2,12 @@
 include("Quadratures.jl")
 @reexport using Gridap.Quadratures
 
+include("DuffyQuadratures.jl")
+@reexport using Gridap.DuffyQuadratures
+
+include("QuadratureFactories.jl")
+@reexport using Gridap.QuadratureFactories
+
 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/QuadratureFactoriesTests.jl

@@ -0,0 +1,17 @@
+module QuadratureFactoriesTests
+
+using Test
+using Gridap
+
+order = 5
+
+p = Polytope(HEX_AXIS,HEX_AXIS)
+quad = Quadrature(p,order)
+@test isa(quad,TensorProductQuadrature)
+
+
+p = Polytope(TET_AXIS,TET_AXIS,TET_AXIS)
+quad = Quadrature(p,order)
+@test isa(quad,DuffyQuadrature)
+
+end # module

test/IntegrationTests/runtests.jl

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