Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

FineToCoarseField for vector-valued fields #970

Merged
merged 7 commits into from
Jan 11, 2024
Merged

FineToCoarseField for vector-valued fields #970

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
706f28d
Generalised FineToCoarseField for vector-valued fields
JordiManyer Jan 10, 2024
4e6afe7
Updated NEWS.md
JordiManyer Jan 10, 2024
3c13ed8
Made FineToCoarseField type stable
JordiManyer Jan 10, 2024
3c34acd
Minor
JordiManyer Jan 10, 2024
16445e1
Remove debug statements
JordiManyer Jan 11, 2024
d49944c
Small bufix
JordiManyer Jan 11, 2024
b583f4d
Small bugfix
JordiManyer Jan 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions NEWS.md
View file
Open in desktop
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
View file
Open in desktop
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
"""
Loading