Skip to content

Commit cea364f

Browse files
fverdugofverdugo
authored
Merge pull request #896 from gridap/add_jacobi_polynomial_bases
Add jacobi polynomial bases
2 parents 8440d03 + 616aa7f commit cea364f

File tree

6 files changed

+429
-1
lines changed

6 files changed

+429
-1
lines changed

‎NEWS.md

+2
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
99

1010
### Added
1111

12+
- Jacobi polynomial bases. Since PR [#896](https://github.com/gridap/Gridap.jl/pull/896).
13+
1214
### Fixed
1315

1416
- Fixed the method `get_normal_vector` for `AdaptedTriangulation`. The method `get_facet_normal`

‎src/Polynomials/JacobiPolynomialBases.jl

+359
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,359 @@
1+
struct JacobiPolynomial <: Field end
2+
3+
struct JacobiPolynomialBasis{D,T} <: AbstractVector{JacobiPolynomial}
4+
orders::NTuple{D,Int}
5+
terms::Vector{CartesianIndex{D}}
6+
function JacobiPolynomialBasis{D}(
7+
::Type{T}, orders::NTuple{D,Int}, terms::Vector{CartesianIndex{D}}) where {D,T}
8+
new{D,T}(orders,terms)
9+
end
10+
end
11+
12+
@inline Base.size(a::JacobiPolynomialBasis{D,T}) where {D,T} = (length(a.terms)*num_components(T),)
13+
@inline Base.getindex(a::JacobiPolynomialBasis,i::Integer) = JacobiPolynomial()
14+
@inline Base.IndexStyle(::JacobiPolynomialBasis) = IndexLinear()
15+
16+
function JacobiPolynomialBasis{D}(
17+
::Type{T}, orders::NTuple{D,Int}, filter::Function=_q_filter) where {D,T}
18+
19+
terms = _define_terms(filter, orders)
20+
JacobiPolynomialBasis{D}(T,orders,terms)
21+
end
22+
23+
function JacobiPolynomialBasis{D}(
24+
::Type{T}, order::Int, filter::Function=_q_filter) where {D,T}
25+
26+
orders = tfill(order,Val{D}())
27+
JacobiPolynomialBasis{D}(T,orders,filter)
28+
end
29+
30+
# API
31+
32+
function get_exponents(b::JacobiPolynomialBasis)
33+
indexbase = 1
34+
[Tuple(t) .- indexbase for t in b.terms]
35+
end
36+
37+
function get_order(b::JacobiPolynomialBasis)
38+
maximum(b.orders)
39+
end
40+
41+
function get_orders(b::JacobiPolynomialBasis)
42+
b.orders
43+
end
44+
45+
return_type(::JacobiPolynomialBasis{D,T}) where {D,T} = T
46+
47+
# Field implementation
48+
49+
function return_cache(f::JacobiPolynomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
50+
@assert D == length(eltype(x)) "Incorrect number of point components"
51+
np = length(x)
52+
ndof = length(f.terms)*num_components(T)
53+
n = 1 + _maximum(f.orders)
54+
r = CachedArray(zeros(T,(np,ndof)))
55+
v = CachedArray(zeros(T,(ndof,)))
56+
c = CachedArray(zeros(eltype(T),(D,n)))
57+
(r, v, c)
58+
end
59+
60+
function evaluate!(cache,f::JacobiPolynomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
61+
r, v, c = cache
62+
np = length(x)
63+
ndof = length(f.terms)*num_components(T)
64+
n = 1 + _maximum(f.orders)
65+
setsize!(r,(np,ndof))
66+
setsize!(v,(ndof,))
67+
setsize!(c,(D,n))
68+
for i in 1:np
69+
@inbounds xi = x[i]
70+
_evaluate_nd_jp!(v,xi,f.orders,f.terms,c)
71+
for j in 1:ndof
72+
@inbounds r[i,j] = v[j]
73+
end
74+
end
75+
r.array
76+
end
77+
78+
function return_cache(
79+
fg::FieldGradientArray{1,JacobiPolynomialBasis{D,V}},
80+
x::AbstractVector{<:Point}) where {D,V}
81+
82+
f = fg.fa
83+
@assert D == length(eltype(x)) "Incorrect number of point components"
84+
np = length(x)
85+
ndof = length(f.terms)*num_components(V)
86+
xi = testitem(x)
87+
T = gradient_type(V,xi)
88+
n = 1 + _maximum(f.orders)
89+
r = CachedArray(zeros(T,(np,ndof)))
90+
v = CachedArray(zeros(T,(ndof,)))
91+
c = CachedArray(zeros(eltype(T),(D,n)))
92+
g = CachedArray(zeros(eltype(T),(D,n)))
93+
(r, v, c, g)
94+
end
95+
96+
function evaluate!(
97+
cache,
98+
fg::FieldGradientArray{1,JacobiPolynomialBasis{D,T}},
99+
x::AbstractVector{<:Point}) where {D,T}
100+
101+
f = fg.fa
102+
r, v, c, g = cache
103+
np = length(x)
104+
ndof = length(f.terms) * num_components(T)
105+
n = 1 + _maximum(f.orders)
106+
setsize!(r,(np,ndof))
107+
setsize!(v,(ndof,))
108+
setsize!(c,(D,n))
109+
setsize!(g,(D,n))
110+
for i in 1:np
111+
@inbounds xi = x[i]
112+
_gradient_nd_jp!(v,xi,f.orders,f.terms,c,g,T)
113+
for j in 1:ndof
114+
@inbounds r[i,j] = v[j]
115+
end
116+
end
117+
r.array
118+
end
119+
120+
function return_cache(
121+
fg::FieldGradientArray{2,JacobiPolynomialBasis{D,V}},
122+
x::AbstractVector{<:Point}) where {D,V}
123+
124+
f = fg.fa
125+
@assert D == length(eltype(x)) "Incorrect number of point components"
126+
np = length(x)
127+
ndof = length(f.terms)*num_components(V)
128+
xi = testitem(x)
129+
T = gradient_type(gradient_type(V,xi),xi)
130+
n = 1 + _maximum(f.orders)
131+
r = CachedArray(zeros(T,(np,ndof)))
132+
v = CachedArray(zeros(T,(ndof,)))
133+
c = CachedArray(zeros(eltype(T),(D,n)))
134+
g = CachedArray(zeros(eltype(T),(D,n)))
135+
h = CachedArray(zeros(eltype(T),(D,n)))
136+
(r, v, c, g, h)
137+
end
138+
139+
function evaluate!(
140+
cache,
141+
fg::FieldGradientArray{2,JacobiPolynomialBasis{D,T}},
142+
x::AbstractVector{<:Point}) where {D,T}
143+
144+
f = fg.fa
145+
r, v, c, g, h = cache
146+
np = length(x)
147+
ndof = length(f.terms) * num_components(T)
148+
n = 1 + _maximum(f.orders)
149+
setsize!(r,(np,ndof))
150+
setsize!(v,(ndof,))
151+
setsize!(c,(D,n))
152+
setsize!(g,(D,n))
153+
setsize!(h,(D,n))
154+
for i in 1:np
155+
@inbounds xi = x[i]
156+
_hessian_nd_jp!(v,xi,f.orders,f.terms,c,g,h,T)
157+
for j in 1:ndof
158+
@inbounds r[i,j] = v[j]
159+
end
160+
end
161+
r.array
162+
end
163+
164+
# Optimizing evaluation at a single point
165+
166+
function return_cache(f::JacobiPolynomialBasis{D,T},x::Point) where {D,T}
167+
ndof = length(f.terms)*num_components(T)
168+
r = CachedArray(zeros(T,(ndof,)))
169+
xs = [x]
170+
cf = return_cache(f,xs)
171+
r, cf, xs
172+
end
173+
174+
function evaluate!(cache,f::JacobiPolynomialBasis{D,T},x::Point) where {D,T}
175+
r, cf, xs = cache
176+
xs[1] = x
177+
v = evaluate!(cf,f,xs)
178+
ndof = size(v,2)
179+
setsize!(r,(ndof,))
180+
a = r.array
181+
copyto!(a,v)
182+
a
183+
end
184+
185+
function return_cache(
186+
f::FieldGradientArray{N,JacobiPolynomialBasis{D,V}}, x::Point) where {N,D,V}
187+
xs = [x]
188+
cf = return_cache(f,xs)
189+
v = evaluate!(cf,f,xs)
190+
r = CachedArray(zeros(eltype(v),(size(v,2),)))
191+
r, cf, xs
192+
end
193+
194+
function evaluate!(
195+
cache, f::FieldGradientArray{N,JacobiPolynomialBasis{D,V}}, x::Point) where {N,D,V}
196+
r, cf, xs = cache
197+
xs[1] = x
198+
v = evaluate!(cf,f,xs)
199+
ndof = size(v,2)
200+
setsize!(r,(ndof,))
201+
a = r.array
202+
copyto!(a,v)
203+
a
204+
end
205+
206+
# Helpers
207+
208+
function _evaluate_1d_jp!(v::AbstractMatrix{T},x,order,d) where T
209+
n = order + 1
210+
z = one(T)
211+
@inbounds v[d,1] = z
212+
if n > 1
213+
ξ = ( 2*x[d] - 1 )
214+
for i in 2:n
215+
@inbounds v[d,i] = sqrt(2*i-1)*jacobi(ξ,i-1,0,0)
216+
end
217+
end
218+
end
219+
220+
function _gradient_1d_jp!(v::AbstractMatrix{T},x,order,d) where T
221+
n = order + 1
222+
z = zero(T)
223+
@inbounds v[d,1] = z
224+
if n > 1
225+
ξ = ( 2*x[d] - 1 )
226+
for i in 2:n
227+
@inbounds v[d,i] = sqrt(2*i-1)*i*jacobi(ξ,i-2,1,1)
228+
end
229+
end
230+
end
231+
232+
function _hessian_1d_jp!(v::AbstractMatrix{T},x,order,d) where T
233+
n = order + 1
234+
z = zero(T)
235+
@inbounds v[d,1] = z
236+
if n > 1
237+
@inbounds v[d,2] = z
238+
ξ = ( 2*x[d] - 1 )
239+
for i in 3:n
240+
@inbounds v[d,i] = sqrt(2*i-1)*(i*(i+1)/2)*jacobi(ξ,i-3,2,2)
241+
end
242+
end
243+
end
244+
245+
function _evaluate_nd_jp!(
246+
v::AbstractVector{V},
247+
x,
248+
orders,
249+
terms::AbstractVector{CartesianIndex{D}},
250+
c::AbstractMatrix{T}) where {V,T,D}
251+
252+
dim = D
253+
for d in 1:dim
254+
_evaluate_1d_jp!(c,x,orders[d],d)
255+
end
256+
257+
o = one(T)
258+
k = 1
259+
260+
for ci in terms
261+
262+
s = o
263+
for d in 1:dim
264+
@inbounds s *= c[d,ci[d]]
265+
end
266+
267+
k = _set_value!(v,s,k)
268+
269+
end
270+
271+
end
272+
273+
function _gradient_nd_jp!(
274+
v::AbstractVector{G},
275+
x,
276+
orders,
277+
terms::AbstractVector{CartesianIndex{D}},
278+
c::AbstractMatrix{T},
279+
g::AbstractMatrix{T},
280+
::Type{V}) where {G,T,D,V}
281+
282+
dim = D
283+
for d in 1:dim
284+
_evaluate_1d_jp!(c,x,orders[d],d)
285+
_gradient_1d_jp!(g,x,orders[d],d)
286+
end
287+
288+
z = zero(Mutable(VectorValue{D,T}))
289+
o = one(T)
290+
k = 1
291+
292+
for ci in terms
293+
294+
s = z
295+
for i in eachindex(s)
296+
@inbounds s[i] = o
297+
end
298+
for q in 1:dim
299+
for d in 1:dim
300+
if d != q
301+
@inbounds s[q] *= c[d,ci[d]]
302+
else
303+
@inbounds s[q] *= g[d,ci[d]]
304+
end
305+
end
306+
end
307+
308+
k = _set_gradient!(v,s,k,V)
309+
310+
end
311+
312+
end
313+
314+
function _hessian_nd_jp!(
315+
v::AbstractVector{G},
316+
x,
317+
orders,
318+
terms::AbstractVector{CartesianIndex{D}},
319+
c::AbstractMatrix{T},
320+
g::AbstractMatrix{T},
321+
h::AbstractMatrix{T},
322+
::Type{V}) where {G,T,D,V}
323+
324+
dim = D
325+
for d in 1:dim
326+
_evaluate_1d_jp!(c,x,orders[d],d)
327+
_gradient_1d_jp!(g,x,orders[d],d)
328+
_hessian_1d_jp!(h,x,orders[d],d)
329+
end
330+
331+
z = zero(Mutable(TensorValue{D,D,T}))
332+
o = one(T)
333+
k = 1
334+
335+
for ci in terms
336+
337+
s = z
338+
for i in eachindex(s)
339+
@inbounds s[i] = o
340+
end
341+
for r in 1:dim
342+
for q in 1:dim
343+
for d in 1:dim
344+
if d != q && d != r
345+
@inbounds s[r,q] *= c[d,ci[d]]
346+
elseif d == q && d ==r
347+
@inbounds s[r,q] *= h[d,ci[d]]
348+
else
349+
@inbounds s[r,q] *= g[d,ci[d]]
350+
end
351+
end
352+
end
353+
end
354+
355+
k = _set_gradient!(v,s,k,V)
356+
357+
end
358+
359+
end

‎src/Polynomials/Polynomials.jl

+3
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@ export QGradMonomialBasis
2525
export QCurlGradMonomialBasis
2626
export PCurlGradMonomialBasis
2727
export ModalC0Basis
28+
export JacobiPolynomialBasis
2829
export get_exponents
2930

3031
export get_order
@@ -41,4 +42,6 @@ include("PCurlGradMonomialBases.jl")
4142

4243
include("ModalC0Bases.jl")
4344

45+
include("JacobiPolynomialBases.jl")
46+
4447
end # module

0 commit comments

Comments
 (0)