Skip to content

Commit

Permalink
Merge pull request #970 from gridap/remove-p4est-warning
Browse files Browse the repository at this point in the history
FineToCoarseField for vector-valued fields
  • Loading branch information
JordiManyer authored Jan 11, 2024
2 parents 10d0b59 + b583f4d commit 5c70c78
Show file tree
Hide file tree
Showing 2 changed files with 59 additions and 97 deletions.
4 changes: 4 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

- Changed how `allocate_vector` works. Now it only allocates, instead of allocating+initialising to zero. Since PR[#963](https://github.com/gridap/Gridap.jl/pull/963).

### Fixed

- Fixed issue where `FineToCoarseField` would not work for vector-valued fields. Since PR[#970](https://github.com/gridap/Gridap.jl/pull/970).

## [0.17.21] - 2023-12-04

### Added
Expand Down
152 changes: 55 additions & 97 deletions src/Adaptivity/FineToCoarseFields.jl
Original file line number Diff line number Diff line change
@@ -1,83 +1,77 @@

"""
Given a domain and a non-overlapping refined cover, a `FineToCoarseField`
is a `Field` defined in the domain and constructed by a set of fields defined on
the subparts of the covering partition.
The refined cover is represented by a `RefinementRule`.
Given a domain and a non-overlapping refined cover, a `FineToCoarseField`
is a `Field` defined in the domain and constructed by a set of fields defined on
the subparts of the covering partition.
The refined cover is represented by a `RefinementRule`.
"""
struct FineToCoarseField{A<:AbstractArray{<:Field},B<:RefinementRule} <: Field
fine_fields :: A
rrule :: B
function FineToCoarseField(fine_fields::AbstractArray{<:Field},rrule::RefinementRule)
is_zero :: Vector{Bool}
function FineToCoarseField(
fine_fields::AbstractArray{<:Field},
rrule::RefinementRule,
is_zero::Vector{Bool}=fill(false,num_subcells(rrule))
)
@check length(fine_fields) == num_subcells(rrule)
A = typeof(fine_fields)
B = typeof(rrule)
new{A,B}(fine_fields,rrule)
new{A,B}(fine_fields,rrule,is_zero)
end
end

# Necessary for distributed meshes, where not all children of a coarse cell may belong to the processor.
function FineToCoarseField(fine_fields::AbstractArray{<:Field},rrule::RefinementRule,child_ids::AbstractArray{<:Integer})
fields = Vector{Field}(undef,num_subcells(rrule))
fields = fill!(fields,ConstantField(0.0))
function FineToCoarseField(
fine_fields::AbstractArray{T},
rrule::RefinementRule,
child_ids::AbstractArray{<:Integer}
) where T <: Field
fields = Vector{Union{T,ZeroField{T}}}(undef,num_subcells(rrule))
fields = fill!(fields,ZeroField(testitem(fine_fields)))
is_zero = fill(true,num_subcells(rrule))
for (k,id) in enumerate(child_ids)
fields[id] = fine_fields[k]
is_zero[id] = false
end
return FineToCoarseField(fields,rrule)
return FineToCoarseField(fields,rrule,is_zero)
end

"""
function Geometry.return_cache(a::FineToCoarseField,x::Point)
fields = a.fine_fields
cmaps = get_inverse_cell_map(a.rrule)
fi_cache = array_cache(fields)
cm_cache = array_cache(cmaps)
xi_cache = Fields.return_cache(first(cmaps),x)
yi_cache = Fields.return_cache(first(fields),x)
return fi_cache, cm_cache, xi_cache, yi_cache
end
function Geometry.evaluate!(cache,a::FineToCoarseField,x::Point)
fi_cache, cm_cache, xi_cache, yi_cache = cache
fields, rr = a.fine_fields, a.rrule
cmaps = get_inverse_cell_map(rr)
child_id = x_to_cell(rr,x) # Find correct subcell
fi = getindex!(fi_cache,fields,child_id)
mi = getindex!(cm_cache,cmaps,child_id)
xi = Fields.evaluate!(xi_cache,mi,x) # xc -> xf
yi = Fields.evaluate!(yi_cache,fi,xi) # xf -> yf
return yi
end
"""

function Geometry.return_cache(a::FineToCoarseField,x::AbstractArray{<:Point})
fields, rr = a.fine_fields, a.rrule
fields, rr, is_zero = a.fine_fields, a.rrule, a.is_zero
cmaps = get_inverse_cell_map(a.rrule)

xi_cache = array_cache(x)
fi_cache = array_cache(fields)
mi_cache = array_cache(cmaps)

xi = getindex!(xi_cache,x,1)
child_id = x_to_cell(rr,xi)
mi = getindex!(mi_cache,cmaps,child_id)
fi = getindex!(fi_cache,fields,child_id)
# Generic caches
child_ids = map(i -> x_to_cell(rr,getindex!(xi_cache,x,i)),eachindex(x))
pos = findfirst(id -> !is_zero[id],child_ids)
xi = getindex!(xi_cache,x,pos)
id = child_ids[pos]
id = x_to_cell(rr,xi)
mi = getindex!(mi_cache,cmaps,id)
fi = getindex!(fi_cache,fields,id)

zi_cache = Fields.return_cache(mi,xi)
zi = evaluate!(zi_cache,mi,xi)

yi_type = Fields.return_type(fi,zi)
yi_cache = Fields.return_cache(fi,zi)
y_cache = Arrays.CachedArray(yi_type,1)
y_cache = Arrays.CachedArray(zeros(yi_type,size(x)))

# Evaluation caches
fi_zero = ZeroField(fi)
yi_nonzero_cache = Fields.return_cache(fi,zi)
yi_zero_cache = Fields.return_cache(fi_zero,zi)
yi_cache = (yi_nonzero_cache,yi_zero_cache)

return fi_cache, mi_cache, xi_cache, zi_cache, yi_cache, y_cache
end

function Geometry.evaluate!(cache,a::FineToCoarseField,x::AbstractArray{<:Point})
fi_cache, mi_cache, xi_cache, zi_cache, yi_cache, y_cache = cache
fields, rr = a.fine_fields, a.rrule
fields, rr, is_zero = a.fine_fields, a.rrule, a.is_zero
cmaps = get_inverse_cell_map(rr)

Arrays.setsize!(y_cache, size(x))
Expand All @@ -88,7 +82,8 @@ function Geometry.evaluate!(cache,a::FineToCoarseField,x::AbstractArray{<:Point}
fi = getindex!(fi_cache,fields,child_id)
mi = getindex!(mi_cache,cmaps,child_id)
zi = Fields.evaluate!(zi_cache,mi,xi)
y_cache.array[i] = Fields.evaluate!(yi_cache,fi,zi)
_yi_cache = yi_cache[is_zero[child_id]+1]
y_cache.array[i] = Fields.evaluate!(_yi_cache,fi,zi)
end
return y_cache.array
end
Expand All @@ -97,33 +92,40 @@ end
# Points are pre-classified into the children cells, which allows for the search to be
# skipped entirely.
function Geometry.return_cache(a::FineToCoarseField,x::AbstractArray{<:Point},child_ids::AbstractArray{<:Integer})
fields = a.fine_fields
fields, is_zero = a.fine_fields, a.is_zero
cmaps = get_inverse_cell_map(a.rrule)

# Generic caches
xi_cache = array_cache(x)
fi_cache = array_cache(fields)
mi_cache = array_cache(cmaps)
id_cache = array_cache(child_ids)

xi = getindex!(xi_cache,x,1)
id = getindex!(id_cache,child_ids,1)
pos = findfirst(id -> !is_zero[id],child_ids)
xi = getindex!(xi_cache,x,pos)
id = getindex!(id_cache,child_ids,pos)
mi = getindex!(mi_cache,cmaps,id)
fi = getindex!(fi_cache,fields,id)

zi_cache = Fields.return_cache(mi,xi)
zi = evaluate!(zi_cache,mi,xi)

yi_type = Fields.return_type(fi,zi)
yi_cache = Fields.return_cache(fi,zi)
y_cache = Arrays.CachedArray(yi_type,1)
y_cache = Arrays.CachedArray(zeros(yi_type,size(x)))

# Evaluation caches
fi_zero = ZeroField(fi)
yi_nonzero_cache = Fields.return_cache(fi,zi)
yi_zero_cache = Fields.return_cache(fi_zero,zi)
yi_cache = (yi_nonzero_cache,yi_zero_cache)

return fi_cache, mi_cache, xi_cache, id_cache, zi_cache, yi_cache, y_cache
end

function Geometry.evaluate!(cache,a::FineToCoarseField,x::AbstractArray{<:Point},child_ids::AbstractArray{<:Integer})
fi_cache, mi_cache, xi_cache, id_cache, zi_cache, yi_cache, y_cache = cache
cmaps = get_inverse_cell_map(a.rrule)
fields = a.fine_fields
fields, is_zero = a.fine_fields, a.is_zero

Arrays.setsize!(y_cache, size(x))

Expand All @@ -133,52 +135,8 @@ function Geometry.evaluate!(cache,a::FineToCoarseField,x::AbstractArray{<:Point}
fi = getindex!(fi_cache,fields,id)
mi = getindex!(mi_cache,cmaps,id)
zi = Fields.evaluate!(zi_cache,mi,xi)
y_cache.array[i] = Fields.evaluate!(yi_cache,fi,zi)
end
return y_cache.array
end

"""
function Geometry.return_cache(a::FineToCoarseField,x::Table{<:Point})
@check length(x) == num_subcells(a.rrule)
fields, rr = a.fine_fields, a.rrule
cmaps = get_inverse_cell_map(a.rrule)
xi_cache = array_cache(x.data)
fi_cache = array_cache(fields)
mi_cache = array_cache(cmaps)
xi = getindex!(xi_cache,x.data,1)
child_id = x_to_cell(rr,xi)
mi = getindex!(mi_cache,cmaps,child_id)
fi = getindex!(fi_cache,fields,child_id)
zi_cache = Fields.return_cache(mi,xi)
zi = evaluate!(zi_cache,mi,xi)
yi_type = Fields.return_type(fi,zi)
yi_cache = Fields.return_cache(fi,zi)
y_cache = Arrays.CachedArray(yi_type,1)
return fi_cache, mi_cache, xi_cache, zi_cache, yi_cache, y_cache
end
function Geometry.evaluate!(cache,a::FineToCoarseField,x::Table{<:Point})
@check length(x) == num_subcells(a.rrule)
fi_cache, mi_cache, xi_cache, zi_cache, yi_cache, y_cache = cache
cmaps = get_inverse_cell_map(a.rrule)
fields = a.fine_fields
Arrays.setsize!(y_cache, size(x.data))
for i in eachindex(x)
for j in x.ptrs[i]:x.ptrs[i+1]
xi = getindex!(xi_cache,x.data,j)
fi = getindex!(fi_cache,fields,i)
mi = getindex!(mi_cache,cmaps,i)
zi = Fields.evaluate!(zi_cache,mi,xi)
y_cache.array[j] = Fields.evaluate!(yi_cache,fi,zi)
end
_yi_cache = yi_cache[is_zero[id]+1]
y_cache.array[i] = Fields.evaluate!(_yi_cache,fi,zi)
end
return y_cache.array
end
"""

0 comments on commit 5c70c78

Please sign in to comment.