Skip to content

Commit 78b27c4

Browse files
fverdugofverdugo
authored
Merge pull request #804 from gridap/nedelec_ho_tets
High order Nedelec elements on simplices
2 parents 26a1886 + 42672a0 commit 78b27c4

File tree

4 files changed

+299
-49
lines changed

4 files changed

+299
-49
lines changed

‎src/Polynomials/QGradMonomialBases.jl

+148-39
Original file line numberDiff line numberDiff line change
@@ -239,22 +239,17 @@ function _gradient_nd_qgrad!(
239239

240240
end
241241

242-
243242
struct NedelecPrebasisOnSimplex{D} <: AbstractVector{Monomial}
244243
order::Int
245244
function NedelecPrebasisOnSimplex{D}(order::Integer) where D
246245
new{D}(Int(order))
247246
end
248247
end
249248

250-
function Base.size(a::NedelecPrebasisOnSimplex{3})
251-
@notimplementedif a.order != 0
252-
(6,)
253-
end
254-
255-
function Base.size(a::NedelecPrebasisOnSimplex{2})
256-
@notimplementedif a.order != 0
257-
(3,)
249+
function Base.size(a::NedelecPrebasisOnSimplex{d}) where d
250+
k = a.order+1
251+
n = div(k*prod(i->(k+i),2:d),factorial(d-1))
252+
(n,)
258253
end
259254

260255
Base.getindex(a::NedelecPrebasisOnSimplex,i::Integer) = Monomial()
@@ -264,52 +259,91 @@ num_terms(a::NedelecPrebasisOnSimplex) = length(a)
264259
get_order(f::NedelecPrebasisOnSimplex) = f.order
265260

266261
function return_cache(
267-
f::NedelecPrebasisOnSimplex{D},x::AbstractVector{<:Point}) where D
268-
@notimplementedif f.order != 0
262+
f::NedelecPrebasisOnSimplex{d},x::AbstractVector{<:Point}) where d
269263
np = length(x)
270264
ndofs = num_terms(f)
271265
V = eltype(x)
272266
a = zeros(V,(np,ndofs))
273-
CachedArray(a)
267+
k = f.order+1
268+
P = MonomialBasis{d}(VectorValue{d,Float64},k-1,(e,order)->sum(e)<=order)
269+
cP = return_cache(P,x)
270+
CachedArray(a), cP, P
274271
end
275272

276273
function evaluate!(
277274
cache,f::NedelecPrebasisOnSimplex{3},x::AbstractVector{<:Point})
278-
@notimplementedif f.order != 0
275+
ca,cP,P = cache
276+
k = f.order+1
279277
np = length(x)
280278
ndofs = num_terms(f)
281-
setsize!(cache,(np,ndofs))
282-
a = cache.array
279+
ndofsP = length(P)
280+
setsize!(ca,(np,ndofs))
281+
Px = evaluate!(cP,P,x)
282+
a = ca.array
283283
V = eltype(x)
284284
T = eltype(V)
285285
z = zero(T)
286286
u = one(T)
287287
for (ip,p) in enumerate(x)
288-
a[ip,1] = VectorValue((u,z,z))
289-
a[ip,2] = VectorValue((z,u,z))
290-
a[ip,3] = VectorValue((z,z,u))
291-
a[ip,4] = VectorValue((-p[2],p[1],z))
292-
a[ip,5] = VectorValue((-p[3],z,p[1]))
293-
a[ip,6] = VectorValue((z,-p[3],p[2]))
288+
for j in 1:ndofsP
289+
a[ip,j] = Px[ip,j]
290+
end
291+
i = ndofsP
292+
x1,x2,x3 = x[ip]
293+
zp = zero(x1)
294+
for β in 1:k
295+
for α in 1:(k+1-β)
296+
i += 1
297+
a[ip,i] = VectorValue(
298+
-x1^(α-1)*x2^(k-α-β+2)*x3^(β-1),
299+
x1^α*x2^(k-α-β+1)*x3^(β-1),
300+
zp)
301+
i += 1
302+
a[ip,i] = VectorValue(
303+
-x1^(k-α-β+1)*x2^(β-1)*x3^α,
304+
zp,
305+
x1^(k-α-β+2)*x2^(β-1)*x3^(α-1))
306+
end
307+
end
308+
for γ in 1:k
309+
i += 1
310+
a[ip,i] = VectorValue(
311+
zp,
312+
-x2^(γ-1)*x3^(k-γ+1),
313+
x2^γ*x3^(k-γ))
314+
end
294315
end
295316
a
296317
end
297318

298319
function evaluate!(
299320
cache,f::NedelecPrebasisOnSimplex{2},x::AbstractVector{<:Point})
300-
@notimplementedif f.order != 0
321+
ca,cP,P = cache
322+
k = f.order+1
301323
np = length(x)
302324
ndofs = num_terms(f)
303-
setsize!(cache,(np,ndofs))
304-
a = cache.array
325+
ndofsP = length(P)
326+
setsize!(ca,(np,ndofs))
327+
a = ca.array
305328
V = eltype(x)
306329
T = eltype(V)
307330
z = zero(T)
308331
u = one(T)
332+
Px = evaluate!(cP,P,x)
309333
for (ip,p) in enumerate(x)
310-
a[ip,1] = VectorValue((u,z))
311-
a[ip,2] = VectorValue((z,u))
312-
a[ip,3] = VectorValue((-p[2],p[1]))
334+
for j in 1:ndofsP
335+
a[ip,j] = Px[ip,j]
336+
end
337+
i = ndofsP
338+
x1,x2 = x[ip]
339+
zp = zero(x1)
340+
for α in 1:k
341+
i += 1
342+
a[ip,i] = VectorValue(-x1^(α-1)*x2^(k-α+1),x1^α*x2^(k-α))
343+
end
344+
#a[ip,1] = VectorValue((u,z))
345+
#a[ip,2] = VectorValue((z,u))
346+
#a[ip,3] = VectorValue((-p[2],p[1]))
313347
end
314348
a
315349
end
@@ -318,56 +352,131 @@ function return_cache(
318352
g::FieldGradientArray{1,<:NedelecPrebasisOnSimplex{D}},
319353
x::AbstractVector{<:Point}) where D
320354
f = g.fa
321-
@notimplementedif f.order != 0
322355
np = length(x)
323356
ndofs = num_terms(f)
324357
xi = testitem(x)
325358
V = eltype(x)
326359
G = gradient_type(V,xi)
327360
a = zeros(G,(np,ndofs))
328-
CachedArray(a)
361+
k = f.order+1
362+
mb = MonomialBasis{D}(VectorValue{D,Float64},k-1,(e,order)->sum(e)<=order)
363+
P = Broadcasting(∇)(mb)
364+
cP = return_cache(P,x)
365+
CachedArray(a), cP, P
329366
end
330367

331368
function evaluate!(
332369
cache,
333370
g::FieldGradientArray{1,<:NedelecPrebasisOnSimplex{3}},
334371
x::AbstractVector{<:Point})
372+
ca,cP,P = cache
335373
f = g.fa
336-
@notimplementedif f.order != 0
374+
k = f.order+1
337375
np = length(x)
338376
ndofs = num_terms(f)
339-
setsize!(cache,(np,ndofs))
340-
a = cache.array
377+
setsize!(ca,(np,ndofs))
378+
a = ca.array
341379
fill!(a,zero(eltype(a)))
380+
ndofsP = length(P)
381+
Px = evaluate!(cP,P,x)
342382
V = eltype(x)
343383
T = eltype(V)
344384
z = zero(T)
345385
u = one(T)
346386
for (ip,p) in enumerate(x)
347-
a[ip,4] = TensorValue((z,-u,z, u,z,z, z,z,z))
348-
a[ip,5] = TensorValue((z,z,-u, z,z,z, u,z,z))
349-
a[ip,6] = TensorValue((z,z,z, z,z,-u, z,u,z))
387+
for j in 1:ndofsP
388+
a[ip,j] = Px[ip,j]
389+
end
390+
i = ndofsP
391+
x1,x2,x3 = x[ip]
392+
zp = zero(x1)
393+
for β in 1:k
394+
for α in 1:(k+1-β)
395+
i += 1
396+
a[ip,i] = TensorValue(
397+
#-x1^(α-1)*x2^(k-α-β+2)*x3^(β-1),
398+
-(α-1)*_exp(x1,α-2)*x2^(k-α-β+2)*x3^(β-1),
399+
-x1^(α-1)*(k-α-β+2)*_exp(x2,k-α-β+1)*x3^(β-1),
400+
-x1^(α-1)*x2^(k-α-β+2)*(β-1)*_exp(x3,β-2),
401+
#x1^α*x2^(k-α-β+1)*x3^(β-1),
402+
α*_exp(x1,α-1)*x2^(k-α-β+1)*x3^(β-1),
403+
x1^α*(k-α-β+1)*_exp(x2,k-α-β)*x3^(β-1),
404+
x1^α*x2^(k-α-β+1)*(β-1)*_exp(x3,β-2),
405+
#zp,
406+
zp,zp,zp)
407+
i += 1
408+
a[ip,i] = TensorValue(
409+
#-x1^(k-α-β+1)*x2^(β-1)*x3^α,
410+
-(k-α-β+1)*_exp(x1,k-α-β)*x2^(β-1)*x3^α,
411+
-x1^(k-α-β+1)*(β-1)*_exp(x2,β-2)*x3^α,
412+
-x1^(k-α-β+1)*x2^(β-1)*α*_exp(x3,α-1),
413+
# zp
414+
zp,zp,zp,
415+
#x1^(k-α-β+2)*x2^(β-1)*x3^(α-1),
416+
(k-α-β+2)*_exp(x1,k-α-β+1)*x2^(β-1)*x3^(α-1),
417+
x1^(k-α-β+2)*(β-1)*_exp(x2,β-2)*x3^(α-1),
418+
x1^(k-α-β+2)*x2^(β-1)*(α-1)*_exp(x3,α-2))
419+
end
420+
end
421+
for γ in 1:k
422+
i += 1
423+
a[ip,i] = TensorValue(
424+
#zp
425+
zp,zp,zp,
426+
#-x2^(γ-1)*x3^(k-γ+1),
427+
-0*x2^(γ-1)*x3^(k-γ+1),
428+
-(γ-1)*_exp(x2,γ-2)*x3^(k-γ+1),
429+
-x2^(γ-1)*(k-γ+1)*_exp(x3,k-γ),
430+
#x2^γ*x3^(k-γ),
431+
0*x2^γ*x3^(k-γ),
432+
γ*_exp(x2,γ-1)*x3^(k-γ),
433+
x2^γ*(k-γ)*_exp(x3,k-γ-1))
434+
end
435+
#a[ip,4] = TensorValue((z,-u,z, u,z,z, z,z,z))
436+
#a[ip,5] = TensorValue((z,z,-u, z,z,z, u,z,z))
437+
#a[ip,6] = TensorValue((z,z,z, z,z,-u, z,u,z))
350438
end
351439
a
352440
end
353441

442+
_exp(a,y) = y>0 ? a^y : one(a)
443+
354444
function evaluate!(
355445
cache,
356446
g::FieldGradientArray{1,<:NedelecPrebasisOnSimplex{2}},
357447
x::AbstractVector{<:Point})
358448
f = g.fa
359-
@notimplementedif f.order != 0
449+
ca,cP,P = cache
450+
k = f.order+1
360451
np = length(x)
361452
ndofs = num_terms(f)
362-
setsize!(cache,(np,ndofs))
363-
a = cache.array
453+
setsize!(ca,(np,ndofs))
454+
a = ca.array
364455
fill!(a,zero(eltype(a)))
365456
V = eltype(x)
366457
T = eltype(V)
367458
z = zero(T)
368459
u = one(T)
460+
ndofsP = length(P)
461+
Px = evaluate!(cP,P,x)
369462
for (ip,p) in enumerate(x)
370-
a[ip,3] = TensorValue((z,-u, u,z))
463+
for j in 1:ndofsP
464+
a[ip,j] = Px[ip,j]
465+
end
466+
i = ndofsP
467+
x1,x2 = x[ip]
468+
zp = zero(x1)
469+
for α in 1:k
470+
i += 1
471+
a[ip,i] = TensorValue(
472+
#-x1^(α-1)*x2^(k-α+1),
473+
-(α-1)*_exp(x1,α-2)*x2^(k-α+1),
474+
-x1^(α-1)*(k-α+1)*_exp(x2,k-α),
475+
#x1^α*x2^(k-α),
476+
α*_exp(x1,α-1)*x2^(k-α),
477+
x1^α*(k-α)*_exp(x2,k-α-1))
478+
end
479+
#a[ip,3] = TensorValue((z,-u, u,z))
371480
end
372481
a
373482
end

‎src/ReferenceFEs/NedelecRefFEs.jl

+68-5
Original file line numberDiff line numberDiff line change
@@ -76,7 +76,7 @@ end
7676

7777
function _Nedelec_nodes_and_moments(::Type{et}, p::Polytope, order::Integer) where et
7878

79-
@notimplementedif !( is_n_cube(p) || (is_simplex(p) && order==0) )
79+
@notimplementedif !( is_n_cube(p) || (is_simplex(p) ) )
8080

8181
D = num_dims(p)
8282
ft = VectorValue{D,et}
@@ -92,14 +92,19 @@ function _Nedelec_nodes_and_moments(::Type{et}, p::Polytope, order::Integer) whe
9292

9393
if ( num_dims(p) == 3 && order > 0)
9494

95-
fcips, fmoments = _Nedelec_face_values(p,et,order)
95+
if is_n_cube(p)
96+
fcips, fmoments = _Nedelec_face_values(p,et,order)
97+
else
98+
fcips, fmoments = _Nedelec_face_values_simplex(p,et,order)
99+
end
100+
96101
frange = get_dimrange(p,D-1)
97102
nf_nodes[frange] = fcips
98103
nf_moments[frange] = fmoments
99104

100105
end
101106

102-
if (order > 0)
107+
if ( is_n_cube(p) && order > 0) || ( is_simplex(p) && order > 1)
103108

104109
ccips, cmoments = _Nedelec_cell_values(p,et,order)
105110
crange = get_dimrange(p,D)
@@ -151,7 +156,7 @@ end
151156
function _Nedelec_face_values(p,et,order)
152157

153158
# Reference facet
154-
@assert is_simplex(p) || is_n_cube(p) "We are assuming that all n-faces of the same n-dim are the same."
159+
@assert is_n_cube(p) "We are assuming that all n-faces of the same n-dim are the same."
155160
fp = Polytope{num_dims(p)-1}(p,1)
156161

157162
# geomap from ref face to polytope faces
@@ -194,6 +199,59 @@ function _Nedelec_face_moments(p, fshfs, c_fips, fcips, fwips)
194199
return cvals
195200
end
196201

202+
function _Nedelec_face_values_simplex(p,et,order)
203+
204+
# Reference facet
205+
@assert is_simplex(p) "We are assuming that all n-faces of the same n-dim are the same."
206+
fp = Polytope{num_dims(p)-1}(p,1)
207+
208+
# geomap from ref face to polytope faces
209+
fgeomap = _ref_face_to_faces_geomap(p,fp)
210+
211+
# Compute integration points at all polynomial edges
212+
degree = (order+1)*2
213+
fquad = Quadrature(fp,degree)
214+
fips = get_coordinates(fquad)
215+
wips = get_weights(fquad)
216+
217+
c_fips, fcips, fwips, fJtips = _nfaces_evaluation_points_weights_with_jac(p, fgeomap, fips, wips)
218+
219+
Df = num_dims(fp)
220+
fshfs = MonomialBasis{Df}(VectorValue{Df,et},order-1,(e,k)->sum(e)<=k)
221+
222+
fmoments = _Nedelec_face_moments_simplex(p, fshfs, c_fips, fcips, fwips, fJtips)
223+
224+
return fcips, fmoments
225+
226+
end
227+
228+
function _nfaces_evaluation_points_weights_with_jac(p, fgeomap, fips, wips)
229+
nc = length(fgeomap)
230+
c_fips = fill(fips,nc)
231+
c_wips = fill(wips,nc)
232+
pquad = lazy_map(evaluate,fgeomap,c_fips)
233+
## Must account for diagonals in simplex discretizations to get the correct
234+
## scaling
235+
Jt1 = lazy_map(∇,fgeomap)
236+
Jt1_ips = lazy_map(evaluate,Jt1,c_fips)
237+
#det_J = lazy_map(Broadcasting(meas),Jt1_ips)
238+
#c_detwips = collect(lazy_map(Broadcasting(*),c_wips,det_J))
239+
c_detwips = c_wips
240+
c_fips, pquad, c_detwips, Jt1_ips
241+
end
242+
243+
function _Nedelec_face_moments_simplex(p, fshfs, c_fips, fcips, fwips, fJtips)
244+
nc = length(c_fips)
245+
cfshfs = fill(fshfs, nc)
246+
cfshfs_fips = lazy_map(evaluate,cfshfs,c_fips)
247+
function weigth(qij,Jti,wi)
248+
Ji = transpose(Jti)
249+
Ji⋅qij*wi
250+
end
251+
cvals = map(Broadcasting(weigth),cfshfs_fips,fJtips,fwips)
252+
return cvals
253+
end
254+
197255
# It provides for every cell the nodes and the moments arrays
198256
function _Nedelec_cell_values(p,et,order)
199257

@@ -204,7 +262,12 @@ function _Nedelec_cell_values(p,et,order)
204262
cwips = get_weights(iquad)
205263

206264
# Cell moments, i.e., M(C)_{ab} = q_C^a(xgp_C^b) w_C^b ⋅ ()
207-
cbasis = QCurlGradMonomialBasis{num_dims(p)}(et,order-1)
265+
if is_n_cube(p)
266+
cbasis = QCurlGradMonomialBasis{num_dims(p)}(et,order-1)
267+
else
268+
D = num_dims(p)
269+
cbasis = MonomialBasis{D}(VectorValue{D,et},order-2,(e,k)->sum(e)<=k)
270+
end
208271
cmoments = _Nedelec_cell_moments(p, cbasis, ccips, cwips )
209272

210273
return [ccips], [cmoments]

0 commit comments

Comments
 (0)