Skip to content

Commit 534b37d

Browse files
fverdugofverdugo
authored
Merge pull request #329 from gridap/better_handling_of_0_length_applied_arrays
Better handling of 0 length applied arrays
2 parents 07009ba + 9a34d6d commit 534b37d

12 files changed

+112
-36
lines changed

‎NEWS.md

+1
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
2626
- Removed struct `DiscRefFE`.
2727

2828
### Fixed
29+
- Better handling of FE terms defined on empty triangulations. Since PR [#329](https://github.com/gridap/Gridap.jl/pull/329).
2930
- Replaced `+=` by `add_entry!`. Since PR [#316](https://github.com/gridap/Gridap.jl/pull/316).
3031
- Minor fix to let Vtk.jl support changes in Vtk 1.7.X versus 1.6.X. Since PR [#324](https://github.com/gridap/Gridap.jl/pull/324).
3132

‎src/Arrays/Apply.jl

+3-17
Original file line numberDiff line numberDiff line change
@@ -173,11 +173,7 @@ end
173173
function AppliedArray(g::AbstractArray,f::AbstractArray...)
174174
gi = testitem(g) #Assumes that all kernels return the same type
175175
fi = testitems(f...)
176-
if length(g) == 0
177-
T = kernel_return_type(gi,fi...)
178-
else
179-
T = typeof(apply_kernel_for_cache(gi,fi...))
180-
end
176+
T = typeof(kernel_testitem(gi,fi...))
181177
AppliedArray(T,g,f...)
182178
end
183179

@@ -219,29 +215,19 @@ function _array_cache(hash,a::AppliedArray)
219215
end
220216
cf = array_caches(hash,a.f...)
221217
cgi = kernel_cache(gi,fi...)
222-
ai = apply_kernel_for_cache!(cgi,gi,fi...)
218+
ai = kernel_testitem!(cgi,gi,fi...)
223219
i = -testitem(eachindex(a))
224220
e = Evaluation((i,),ai)
225221
c = (cg, cgi, cf)
226222
(c,e)
227223
end
228224

229-
function apply_kernel_for_cache(k,x...)
230-
cache = kernel_cache(k,x...)
231-
apply_kernel_for_cache!(cache,k,x...)
232-
end
233-
234-
@inline function apply_kernel_for_cache!(cache,k,x...)
235-
apply_kernel!(cache,k,x...)
236-
end
237-
238225
function testitem(a::AppliedArray)
239226
cg = array_cache(a.g)
240227
gi = testitem(a.g)
241228
fi = testitems(a.f...)
242229
cf = array_caches(a.f...)
243-
cgi = kernel_cache(gi,fi...)
244-
ai = apply_kernel_for_cache!(cgi,gi,fi...)
230+
ai = kernel_testitem(gi,fi...)
245231
ai
246232
end
247233

‎src/Arrays/Arrays.jl

+1
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,7 @@ export elem
5252
export contract
5353
export kernel_return_type
5454
export kernel_return_types
55+
export kernel_testitem
5556
export Kernel
5657

5758
export apply

‎src/Arrays/Kernels.jl

+11
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,8 @@ function test_kernel(f,x::Tuple,y,cmp=(==))
8080
@test cmp(z,y)
8181
z = apply_kernel!(cache,f,x...)
8282
@test cmp(z,y)
83+
z = kernel_testitem!(cache,f,x...)
84+
@test cmp(typeof(z),typeof(y))
8385
end
8486

8587

@@ -173,6 +175,15 @@ function _kernel_return_types(x::Tuple,a)
173175
(Ta,)
174176
end
175177

178+
function kernel_testitem(k,x...)
179+
cache = kernel_cache(k,x...)
180+
kernel_testitem!(cache,k,x...)
181+
end
182+
183+
@inline function kernel_testitem!(cache,k,x...)
184+
apply_kernel!(cache,k,x...)
185+
end
186+
176187
# Include some well-known types in this interface
177188

178189
function kernel_return_type(f::Function,x...)

‎src/FESpaces/FESpaces.jl

+1-1
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@ import Gridap.Arrays: getindex!
3838
import Gridap.Arrays: kernel_cache
3939
import Gridap.Arrays: apply_kernel!
4040
import Gridap.Arrays: kernel_return_type
41-
import Gridap.Arrays: apply_kernel_for_cache!
41+
import Gridap.Arrays: kernel_testitem!
4242
import Gridap.Arrays: reindex
4343
import Gridap.Arrays: apply
4444
import Gridap.Geometry: get_cell_map

‎src/FESpaces/StateLaws.jl

+1-1
Original file line numberDiff line numberDiff line change
@@ -63,7 +63,7 @@ function apply_kernel!(cache,k::StateLawKernel,a::AbstractVector,b::AbstractVect
6363
v
6464
end
6565

66-
function apply_kernel_for_cache!(cache,k::StateLawKernel,a::AbstractVector,b::AbstractVector...)
66+
function kernel_testitem!(cache,k::StateLawKernel,a::AbstractVector,b::AbstractVector...)
6767
Q = length(a)
6868
setsize!(cache,size(a))
6969
v = cache.array

‎src/Fields/FieldOperations.jl

+56-8
Original file line numberDiff line numberDiff line change
@@ -82,7 +82,6 @@ end
8282

8383
function kernel_cache(
8484
k::FieldBinOp,a::AbstractVector,b::AbstractVector)
85-
_field_bin_op_checks_vecvec(a,b)
8685
na = length(a)
8786
Ta = eltype(a)
8887
Tb = eltype(b)
@@ -91,6 +90,15 @@ function kernel_cache(
9190
CachedArray(r)
9291
end
9392

93+
function kernel_testitem!(
94+
c,k::FieldBinOp,a::AbstractVector,b::AbstractVector)
95+
if _valid_checks_vecvec(a,b)
96+
apply_kernel!(c,k,a,b)
97+
else
98+
c.array
99+
end
100+
end
101+
94102
@inline function apply_kernel!(
95103
c,k::FieldBinOp,a::AbstractVector,b::AbstractVector)
96104
_field_bin_op_checks_vecvec(a,b)
@@ -104,23 +112,35 @@ end
104112
end
105113

106114
function _field_bin_op_checks_vecvec(a,b)
115+
@assert _valid_checks_vecvec(a,b) "Binary operation between fields: vector vs vector size mismatch."
116+
end
117+
118+
function _valid_checks_vecvec(a,b)
107119
na = length(a)
108120
nb = length(b)
109-
@assert na == nb "Binary operation between fields: vector vs vector size mismatch."
121+
na == nb
110122
end
111123

112124
# Matrix vs Vector
113125

114126
function kernel_cache(
115127
k::FieldBinOp,a::AbstractMatrix,b::AbstractVector)
116-
_field_bin_op_checks_matvec(a,b)
117128
Ta = eltype(a)
118129
Tb = eltype(b)
119130
T = return_type(k.op,Ta,Tb)
120131
r = zeros(T,size(a))
121132
CachedArray(r)
122133
end
123134

135+
function kernel_testitem!(
136+
c,k::FieldBinOp,a::AbstractMatrix,b::AbstractVector)
137+
if _valid_checks_matvec(a,b)
138+
apply_kernel!(c,k,a,b)
139+
else
140+
c.array
141+
end
142+
end
143+
124144
@inline function apply_kernel!(
125145
c,k::FieldBinOp,a::AbstractMatrix,b::AbstractVector)
126146
_field_bin_op_checks_matvec(a,b)
@@ -137,23 +157,35 @@ end
137157
end
138158

139159
function _field_bin_op_checks_matvec(a,b)
160+
@assert _valid_checks_matvec(a,b) "Binary operation between fields: matrix vs vector size mismatch."
161+
end
162+
163+
function _valid_checks_matvec(a,b)
140164
na, _ = size(a)
141165
nb = length(b)
142-
@assert na == nb "Binary operation between fields: matrix vs vector size mismatch."
166+
na == nb
143167
end
144168

145169
# Vector vs Matrix
146170

147171
function kernel_cache(
148172
k::FieldBinOp,a::AbstractVector,b::AbstractMatrix)
149-
_field_bin_op_checks_vecmat(a,b)
150173
Ta = eltype(a)
151174
Tb = eltype(b)
152175
T = return_type(k.op,Ta,Tb)
153176
r = zeros(T,size(b))
154177
CachedArray(r)
155178
end
156179

180+
function kernel_testitem!(
181+
c,k::FieldBinOp,a::AbstractVector,b::AbstractMatrix)
182+
if _valid_checks_vecmat(a,b)
183+
apply_kernel!(c,k,a,b)
184+
else
185+
c.array
186+
end
187+
end
188+
157189
@inline function apply_kernel!(
158190
c,k::FieldBinOp,a::AbstractVector,b::AbstractMatrix)
159191
_field_bin_op_checks_vecmat(a,b)
@@ -170,16 +202,19 @@ end
170202
end
171203

172204
function _field_bin_op_checks_vecmat(a,b)
205+
@assert _valid_checks_vecmat(a,b) "Binary operation between fields: vector vs matrix size mismatch."
206+
end
207+
208+
function _valid_checks_vecmat(a,b)
173209
nb, _ = size(b)
174210
na = length(a)
175-
@assert na == nb "Binary operation between fields: vector vs matrix size mismatch."
211+
na == nb
176212
end
177213

178214
# Matrix vs matrix
179215

180216
function kernel_cache(
181217
k::FieldBinOp,a::AbstractMatrix,b::AbstractMatrix)
182-
_field_bin_op_checks_matmat(a,b)
183218
Ta = eltype(a)
184219
Tb = eltype(b)
185220
T = return_type(k.op,Ta,Tb)
@@ -189,6 +224,15 @@ function kernel_cache(
189224
CachedArray(r)
190225
end
191226

227+
function kernel_testitem!(
228+
c,k::FieldBinOp,a::AbstractMatrix,b::AbstractMatrix)
229+
if _valid_checks_matmat(a,b)
230+
apply_kernel!(c,k,a,b)
231+
else
232+
c.array
233+
end
234+
end
235+
192236
@inline function apply_kernel!(
193237
c,k::FieldBinOp,a::AbstractMatrix,b::AbstractMatrix)
194238
_field_bin_op_checks_matmat(a,b)
@@ -207,9 +251,13 @@ end
207251
end
208252

209253
function _field_bin_op_checks_matmat(a,b)
254+
@assert _valid_checks_matmat(a,b) "Binary operation between fields: matrix vs matrix size mismatch."
255+
end
256+
257+
function _valid_checks_matmat(a,b)
210258
na, ni = size(a)
211259
nb, nj = size(b)
212-
@assert na == nb "Binary operation between fields: matrix vs matrix size mismatch."
260+
na == nb
213261
end
214262

215263
# Define gradients

‎src/Fields/Fields.jl

+1
Original file line numberDiff line numberDiff line change
@@ -76,6 +76,7 @@ export function_field
7676
import Gridap.Arrays: kernel_cache
7777
import Gridap.Arrays: apply_kernel!
7878
import Gridap.Arrays: kernel_return_type
79+
import Gridap.Arrays: kernel_testitem!
7980
import Gridap.TensorValues: outer
8081
import Gridap.TensorValues: inner
8182
import Gridap.TensorValues: symmetric_part

‎src/Fields/Integrate.jl

+21-3
Original file line numberDiff line numberDiff line change
@@ -23,11 +23,18 @@ end
2323
struct IntKernel <: Kernel end
2424

2525
function kernel_cache(k::IntKernel,f::AbstractVector,w,j)
26-
_integrate_checks(f,w,j)
2726
T = _integrate_rt(f,w,j)
2827
zero(T)
2928
end
3029

30+
function kernel_testitem!(c,k::IntKernel,f::AbstractVector,w,j)
31+
if _integrate_valid_sizes(f,w,j)
32+
apply_kernel!(c,k,f,w,j)
33+
else
34+
c
35+
end
36+
end
37+
3138
@noinline function apply_kernel!(z,k::IntKernel,f::AbstractVector,w,j)
3239
_integrate_checks(f,w,j)
3340
r = z
@@ -38,13 +45,20 @@ end
3845
end
3946

4047
function kernel_cache(k::IntKernel,f::AbstractArray,w,j)
41-
_integrate_checks(f,w,j)
4248
T = _integrate_rt(f,w,j)
4349
_, s = _split(size(f)...)
4450
r = zeros(T,s)
4551
c = CachedArray(r)
4652
end
4753

54+
function kernel_testitem!(c,k::IntKernel,f::AbstractArray,w,j)
55+
if _integrate_valid_sizes(f,w,j)
56+
apply_kernel!(c,k,f,w,j)
57+
else
58+
c.array
59+
end
60+
end
61+
4862
@inline function apply_kernel!(c,k::IntKernel,f::AbstractArray,w,j)
4963
_integrate_checks(f,w,j)
5064
np, s = _split(size(f)...)
@@ -72,9 +86,13 @@ function _integrate_rt(f,w,j)
7286
end
7387

7488
function _integrate_checks(f,w,j)
89+
@assert _integrate_valid_sizes(f,w,j) "integrate: sizes mismatch."
90+
end
91+
92+
function _integrate_valid_sizes(f,w,j)
7593
nf, = size(f)
7694
nw = length(w)
7795
nj = length(j)
78-
@assert (nf == nw) && (nw == nj) "integrate: sizes mismatch."
96+
(nf == nw) && (nw == nj)
7997
end
8098

‎src/Fields/MockFields.jl

+3-3
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ function evaluate_field!(c,f::MockField,x)
1919
@inbounds xi = x[i]
2020
@inbounds c[i] = f.v*xi[1]
2121
end
22-
c
22+
c.array
2323
end
2424

2525
@inline function field_gradient(f::MockField{T,D}) where {T,D}
@@ -57,7 +57,7 @@ function evaluate_field!(v,f::MockBasis,x)
5757
@inbounds v[i,j] = f.v*xi[1]
5858
end
5959
end
60-
v
60+
v.array
6161
end
6262

6363
@inline function field_gradient(f::MockBasis{T,D}) where {T,D}
@@ -94,7 +94,7 @@ function evaluate_field!(v,f::OtherMockBasis,x)
9494
@inbounds v[i,j] = 2*xi
9595
end
9696
end
97-
v
97+
v.array
9898
end
9999

100100
@inline function field_gradient(f::OtherMockBasis{D}) where D

‎src/Polynomials/QGradMonomialBases.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -69,7 +69,7 @@ function evaluate_field!(cache,f::QGradMonomialBasis{D,T},x) where {D,T}
6969
@inbounds r[i,j] = v[j]
7070
end
7171
end
72-
r
72+
r.array
7373
end
7474

7575
function gradient_cache(f::QGradMonomialBasis{D,T},x) where {D,T}
@@ -104,7 +104,7 @@ function evaluate_gradient!(cache,f::QGradMonomialBasis{D,T},x) where {D,T}
104104
@inbounds r[i,j] = v[j]
105105
end
106106
end
107-
r
107+
r.array
108108
end
109109

110110
# Helpers

‎test/GridapTests/StokesNitscheTests.jl

+11-1
Original file line numberDiff line numberDiff line change
@@ -103,8 +103,18 @@ t_Ω = AffineFETerm(A_Ω,B_Ω,trian,quad)
103103
# t_∂Ω = FESource(B_∂Ω,btrian,bquad)
104104
t_∂Ω = AffineFETerm(A_∂Ω,B_∂Ω,btrian,bquad)
105105

106+
# Dummy term with 0 facets to check degenerated case
107+
trian0 = BoundaryTriangulation(model,fill(false,Gridap.ReferenceFEs.num_facets(model)))
108+
quad0 = CellQuadrature(trian0,2*order)
109+
s0(x) = VectorValue(2.0*x[1],2.0)
110+
function L0(y)
111+
v,q = y
112+
s0⋅v
113+
end
114+
t_0 = FESource(L0,trian0,quad0)
115+
106116
# Define the FEOperator
107-
op = AffineFEOperator(X,Y,t_Ω,t_∂Ω)
117+
op = AffineFEOperator(X,Y,t_Ω,t_∂Ω,t_0)
108118
# op = LinearFEOperator(Yh,Xh,t_Ω)
109119

110120
# Solve!

0 commit comments

Comments
 (0)