From 8cc891862fd7657b0440a02f7115d392dc10c185 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 30 Oct 2024 00:00:08 +1100 Subject: [PATCH 1/3] AffineMap is now AffineField. Introduced new AffineMap and ConstantMap. --- src/Adaptivity/RefinementRules.jl | 2 +- src/Fields/AffineMaps.jl | 52 ++++++++++++++----- src/Fields/Fields.jl | 2 +- src/Fields/FieldsInterfaces.jl | 10 ++++ src/Geometry/CartesianGrids.jl | 6 +-- test/FieldsTests/AffineMapsTests.jl | 4 +- test/FieldsTests/InverseFieldsTests.jl | 8 +-- .../BoundaryTriangulationsTests.jl | 2 +- test/GridapTests/issue_879.jl | 2 +- 9 files changed, 61 insertions(+), 27 deletions(-) diff --git a/src/Adaptivity/RefinementRules.jl b/src/Adaptivity/RefinementRules.jl index a93d93d20..d617cb091 100644 --- a/src/Adaptivity/RefinementRules.jl +++ b/src/Adaptivity/RefinementRules.jl @@ -25,7 +25,7 @@ end # - The reason why we are saving both the cell maps and the inverse cell maps is to avoid recomputing # them when needed. This is needed for performance when the RefinementRule is used for MacroFEs. # Also, in the case the ref_grid comes from a CartesianGrid, we save the cell maps as -# AffineMaps, which are more efficient than the default linear_combinations. +# AffineFields, which are more efficient than the default linear_combinations. # - We cannot parametrise the RefinementRule by all it's fields, because we will have different types of # RefinementRules in a single mesh. It's the same reason why we don't parametrise the ReferenceFE type. diff --git a/src/Fields/AffineMaps.jl b/src/Fields/AffineMaps.jl index 9ce5725fc..d8f521b3d 100644 --- a/src/Fields/AffineMaps.jl +++ b/src/Fields/AffineMaps.jl @@ -1,12 +1,28 @@ +struct AffineMap <: Map end + +function return_cache(::AffineMap,G::TensorValue{D1,D2,T,L},y0::Point{D2,T},x::Point{D1,T}) + nothing +end + +function evaluate!( + cache,::AffineMap,G::TensorValue{D1,D2},y0::Point{D2},x::Point{D1} +) where {D1,D2} + x⋅G + y0 +end + +struct InverseAffineMap <: Map end + + + """ A Field with this form y = x⋅G + y0 """ -struct AffineMap{D1,D2,T,L} <:Field +struct AffineField{D1,D2,T,L} <: Field gradient::TensorValue{D1,D2,T,L} origin::Point{D2,T} - function AffineMap( + function AffineField( gradient::TensorValue{D1,D2,T,L}, origin::Point{D2,T}) where {D1,D2,T,L} @@ -14,21 +30,21 @@ struct AffineMap{D1,D2,T,L} <:Field end end -affine_map(gradient,origin) = AffineMap(gradient,origin) +affine_map(gradient,origin) = AffineField(gradient,origin) -function evaluate!(cache,f::AffineMap,x::Point) +function evaluate!(cache,f::AffineField,x::Point) G = f.gradient y0 = f.origin x⋅G + y0 end -function return_cache(f::AffineMap,x::AbstractVector{<:Point}) +function return_cache(f::AffineField,x::AbstractVector{<:Point}) T = return_type(f,testitem(x)) y = similar(x,T,size(x)) CachedArray(y) end -function evaluate!(cache,f::AffineMap,x::AbstractVector{<:Point}) +function evaluate!(cache,f::AffineField,x::AbstractVector{<:Point}) setsize!(cache,size(x)) y = cache.array G = f.gradient @@ -41,11 +57,11 @@ function evaluate!(cache,f::AffineMap,x::AbstractVector{<:Point}) y end -function gradient(h::AffineMap) +function gradient(h::AffineField) ConstantField(h.gradient) end -function push_∇∇(∇∇a::Field,ϕ::AffineMap) +function push_∇∇(∇∇a::Field,ϕ::AffineField) # Assuming ϕ is affine map Jt = ∇(ϕ) Jt_inv = pinvJt(Jt) @@ -80,7 +96,7 @@ end function lazy_map( k::Broadcasting{typeof(push_∇∇)}, cell_∇∇a::AbstractArray, - cell_map::AbstractArray{<:AffineMap}) + cell_map::AbstractArray{<:AffineField}) cell_Jt = lazy_map(∇,cell_map) cell_invJt = lazy_map(Operation(pinvJt),cell_Jt) lazy_map(Broadcasting(Operation(push_∇∇)),cell_∇∇a,cell_invJt) @@ -89,19 +105,19 @@ end function lazy_map( k::Broadcasting{typeof(push_∇∇)}, cell_∇∇at::LazyArray{<:Fill{typeof(transpose)}}, - cell_map::AbstractArray{<:AffineMap}) + cell_map::AbstractArray{<:AffineField}) cell_∇∇a = cell_∇∇at.args[1] cell_∇∇b = lazy_map(k,cell_∇∇a,cell_map) cell_∇∇bt = lazy_map(transpose,cell_∇∇b) cell_∇∇bt end -function inverse_map(f::AffineMap) +function inverse_map(f::AffineField) Jt = f.gradient y0 = f.origin invJt = pinvJt(Jt) x0 = -y0⋅invJt - AffineMap(invJt,x0) + AffineField(invJt,x0) end function lazy_map(::typeof(∇),a::LazyArray{<:Fill{typeof(affine_map)}}) @@ -109,8 +125,16 @@ function lazy_map(::typeof(∇),a::LazyArray{<:Fill{typeof(affine_map)}}) lazy_map(constant_field,gradients) end -function Base.zero(::Type{<:AffineMap{D1,D2,T}}) where {D1,D2,T} +function lazy_map( + ::typeof(evaluate),a::LazyArray{<:Fill{typeof(affine_map)}},x::AbstractArray +) + gradients = a.args[1] + origins = a.args[2] + lazy_map(Broadcasting(AffineMap()),gradients,origins,x) +end + +function Base.zero(::Type{<:AffineField{D1,D2,T}}) where {D1,D2,T} gradient = TensorValue{D1,D2}(tfill(zero(T),Val{D1*D2}())) origin = Point{D2,T}(tfill(zero(T),Val{D2}())) - AffineMap(gradient,origin) + AffineField(gradient,origin) end diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl index 05e5e7a56..0dd060b4c 100644 --- a/src/Fields/Fields.jl +++ b/src/Fields/Fields.jl @@ -56,7 +56,7 @@ export MockFieldArray export Point export inverse_map -export AffineMap +export AffineField export affine_map export gradient diff --git a/src/Fields/FieldsInterfaces.jl b/src/Fields/FieldsInterfaces.jl index bdad6a111..799bc96cc 100644 --- a/src/Fields/FieldsInterfaces.jl +++ b/src/Fields/FieldsInterfaces.jl @@ -288,6 +288,16 @@ function lazy_map(::Operation{typeof(inv)},a::LazyArray{<:Fill{typeof(constant_f lazy_map(constant_field,vinv) end +struct ConstantMap <: Map end + +return_cache(::ConstantMap,v::Number,x::Point) = nothing +evaluate!(cache,::ConstantMap,v::Number,x::Point) = v + +function lazy_map(::typeof(evaluate),a::LazyArray{<:Fill{typeof(constant_field)}},x::AbstractArray) + values = a.args[1] + lazy_map(Broadcasting(ConstantMap()),values,x) +end + ## Make Function behave like Field return_cache(f::FieldGradient{N,<:Function},x::Point) where N = gradient(f.object,Val(N)) diff --git a/src/Geometry/CartesianGrids.jl b/src/Geometry/CartesianGrids.jl index ee83694c4..5fd239957 100644 --- a/src/Geometry/CartesianGrids.jl +++ b/src/Geometry/CartesianGrids.jl @@ -237,7 +237,7 @@ end # Cell map -struct CartesianMap{D,T,L} <: AbstractArray{AffineMap{D,D,T,L},D} +struct CartesianMap{D,T,L} <: AbstractArray{AffineField{D,D,T,L},D} data::CartesianDescriptor{D,T,typeof(identity)} function CartesianMap(des::CartesianDescriptor{D,T}) where {D,T} L = D*D @@ -256,9 +256,9 @@ function Base.getindex(a::CartesianMap{D,T},I::Vararg{Integer,D}) where {D,T} @inbounds for d in 1:D p[d] = x0[d] + (I[d]-1)*dx[d] end - origin = Point(p) + origin = Point(p) grad = diagonal_tensor(VectorValue(dx)) - AffineMap(grad,origin) + AffineField(grad,origin) end function lazy_map(::typeof(∇),a::CartesianMap) diff --git a/test/FieldsTests/AffineMapsTests.jl b/test/FieldsTests/AffineMapsTests.jl index 284aaec1d..963e3faf6 100644 --- a/test/FieldsTests/AffineMapsTests.jl +++ b/test/FieldsTests/AffineMapsTests.jl @@ -9,7 +9,7 @@ using Test origin = Point(1,1) g = TensorValue(2,0,0,2) -h = AffineMap(g,origin) +h = AffineField(g,origin) @test isa(∇(h),ConstantField) @test isa(Broadcasting(∇)(h),ConstantField) @@ -39,7 +39,7 @@ cell_to_∇hx = lazy_map(evaluate,cell_to_∇h,cell_to_x) test_array(cell_to_hx,fill(hx,ncells)) test_array(cell_to_∇hx,fill(∇hx,ncells)) -T = AffineMap{3,3,Int} +T = AffineField{3,3,Int} @test isa(zero(T),T) #display(cell_to_hx) diff --git a/test/FieldsTests/InverseFieldsTests.jl b/test/FieldsTests/InverseFieldsTests.jl index 8cf7bb5d2..f8db6750a 100644 --- a/test/FieldsTests/InverseFieldsTests.jl +++ b/test/FieldsTests/InverseFieldsTests.jl @@ -8,22 +8,22 @@ using Test b0 = Point(0,0) m0 = TensorValue(1,0,0,1) -id = AffineMap(m0,b0) +id = AffineField(m0,b0) b1 = Point(1,1) m1 = TensorValue(2,0,0,2) -h1 = AffineMap(m1,b1) +h1 = AffineField(m1,b1) b2 = Point(3,3) m2 = TensorValue(4,0,0,4) -h2 = AffineMap(m2,b2) +h2 = AffineField(m2,b2) h = h2 ∘ h1 # (4x+3) ∘ (2x+1) = 8x+7 b3 = Point(7,7) m3 = TensorValue(8,0,0,8) -h3 = AffineMap(m3,b3) +h3 = AffineField(m3,b3) h1inv = inverse_map(h1) h2inv = inverse_map(h2) diff --git a/test/GeometryTests/BoundaryTriangulationsTests.jl b/test/GeometryTests/BoundaryTriangulationsTests.jl index f876d2a86..b775798ef 100644 --- a/test/GeometryTests/BoundaryTriangulationsTests.jl +++ b/test/GeometryTests/BoundaryTriangulationsTests.jl @@ -78,7 +78,7 @@ glue = get_glue(btrian,Val(1)) glue = get_glue(btrian,Val(2)) @test glue.tface_to_mface === btrian.glue.face_to_cell -@test isa(get_cell_map(btrian)[1],AffineMap) +@test isa(get_cell_map(btrian)[1],AffineField) face_s_q = glue.tface_to_mface_map diff --git a/test/GridapTests/issue_879.jl b/test/GridapTests/issue_879.jl index d0eefea69..bd9613aad 100644 --- a/test/GridapTests/issue_879.jl +++ b/test/GridapTests/issue_879.jl @@ -38,6 +38,6 @@ fΓ = interpolate_everywhere(fa, FESpace(Γ,reffe)) # ERROR: LoadError: DimensionMismatch: matrix is not square: dimensions are (1, 2) # Corrections -# Modified src/Fields/AffineMaps.jl +# Modified src/Fields/AffineFields.jl end \ No newline at end of file From 028897471244ae313d32e649d8d6dfd58f2daee5 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 30 Oct 2024 00:04:33 +1100 Subject: [PATCH 2/3] Minor fix --- src/Fields/AffineMaps.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/Fields/AffineMaps.jl b/src/Fields/AffineMaps.jl index d8f521b3d..e99221107 100644 --- a/src/Fields/AffineMaps.jl +++ b/src/Fields/AffineMaps.jl @@ -1,7 +1,7 @@ struct AffineMap <: Map end -function return_cache(::AffineMap,G::TensorValue{D1,D2,T,L},y0::Point{D2,T},x::Point{D1,T}) +function return_cache(::AffineMap,G::TensorValue{D1,D2},y0::Point{D2},x::Point{D1}) where {D1,D2} nothing end @@ -11,9 +11,6 @@ function evaluate!( x⋅G + y0 end -struct InverseAffineMap <: Map end - - """ A Field with this form From dfd7af5e008446ffb34a709f27d62a6991a0dee4 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 30 Oct 2024 16:31:52 +1100 Subject: [PATCH 3/3] Optimizations for lazy_map(evaluate,lazy_map(Reindex(f),ids),x) --- src/Arrays/Reindex.jl | 8 ++++++++ src/Fields/ApplyOptimizations.jl | 4 ++-- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/Arrays/Reindex.jl b/src/Arrays/Reindex.jl index 7c6442790..429d5b8b3 100644 --- a/src/Arrays/Reindex.jl +++ b/src/Arrays/Reindex.jl @@ -129,3 +129,11 @@ end i = j_to_i[j] i_to_v[i]=v end + +# This optimization is important when integrating on partial domains (Triangulation views) +function lazy_map(::typeof(evaluate),a::LazyArray{<:Fill{<:Reindex}},x::AbstractArray) + fields = a.maps.value.values + ids = a.args[1] + vals = lazy_map(evaluate,fields,x) + return lazy_map(Reindex(vals),ids) +end diff --git a/src/Fields/ApplyOptimizations.jl b/src/Fields/ApplyOptimizations.jl index 76dd091fd..95df7c730 100644 --- a/src/Fields/ApplyOptimizations.jl +++ b/src/Fields/ApplyOptimizations.jl @@ -135,14 +135,14 @@ function lazy_map( k::Broadcasting{typeof(∇)}, a::LazyArray{<:Fill{typeof(transpose)}}) i_to_basis = lazy_map(k,a.args[1]) - lazy_map( transpose, i_to_basis) + lazy_map(transpose, i_to_basis) end function lazy_map( k::Broadcasting{typeof(∇∇)}, a::LazyArray{<:Fill{typeof(transpose)}}) i_to_basis = lazy_map(k,a.args[1]) - lazy_map( transpose, i_to_basis) + lazy_map(transpose, i_to_basis) end # Gradient rules