diff --git a/NEWS.md b/NEWS.md index e7e171a2d..62100f716 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,12 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## Unreleased + +### Fixed + +- `BlockMultiFieldStyle` available for `TransientMultiFieldFESpaces` since PR [#946](https://github.com/gridap/Gridap.jl/pull/946). + ## [0.17.20] - 2023-10-01 ### Added diff --git a/src/MultiField/BlockSparseMatrixAssemblers.jl b/src/MultiField/BlockSparseMatrixAssemblers.jl index 105cb910d..b4a336d74 100644 --- a/src/MultiField/BlockSparseMatrixAssemblers.jl +++ b/src/MultiField/BlockSparseMatrixAssemblers.jl @@ -80,8 +80,9 @@ function BlockSparseMatrixAssembler(trial::MultiFieldFESpace, @notimplemented msg end -function BlockSparseMatrixAssembler(trial::MultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}, - test::MultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}, +function BlockSparseMatrixAssembler(::BlockMultiFieldStyle{NB,SB,P}, + trial, + test, matrix_builder, vector_builder, strategy=FESpaces.DefaultAssemblyStrategy()) where {NB,SB,P} @@ -103,12 +104,13 @@ function BlockSparseMatrixAssembler(trial::MultiFieldFESpace{<:BlockMultiFieldSt return BlockSparseMatrixAssembler{NB,NV,SB,P}(block_assemblers) end -function FESpaces.SparseMatrixAssembler(mat, - vec, - trial::MultiFieldFESpace{<:BlockMultiFieldStyle}, - test ::MultiFieldFESpace{<:BlockMultiFieldStyle}, - strategy::AssemblyStrategy=DefaultAssemblyStrategy()) - return BlockSparseMatrixAssembler(trial,test,SparseMatrixBuilder(mat),ArrayBuilder(vec),strategy) +function FESpaces.SparseMatrixAssembler(mat,vec, + trial::MultiFieldFESpace{MS}, + test ::MultiFieldFESpace{MS}, + strategy::AssemblyStrategy=DefaultAssemblyStrategy() + ) where MS <: BlockMultiFieldStyle + mfs = MultiFieldStyle(test) + return BlockSparseMatrixAssembler(mfs,trial,test,SparseMatrixBuilder(mat),ArrayBuilder(vec),strategy) end # BlockArrays extensions @@ -224,7 +226,7 @@ for T in (:AddEntriesMap,:TouchEntriesMap) cache end - function Fields.evaluate!(cache, k::$T,A::$MT,v::MatrixBlock,I::VectorBlock,J::VectorBlock) + function Fields.evaluate!(cache,k::$T,A::$MT,v::MatrixBlock,I::VectorBlock,J::VectorBlock) ni,nj = size(v.touched) for j in 1:nj for i in 1:ni diff --git a/src/MultiField/MultiFieldFESpaces.jl b/src/MultiField/MultiFieldFESpaces.jl index 436e28979..0df1528a7 100644 --- a/src/MultiField/MultiFieldFESpaces.jl +++ b/src/MultiField/MultiFieldFESpaces.jl @@ -34,12 +34,12 @@ function BlockMultiFieldStyle(NB::Integer) return BlockMultiFieldStyle(NB,SB) end -function BlockMultiFieldStyle(::BlockMultiFieldStyle{NB,SB,P},spaces::Vector{<:SingleFieldFESpace}) where {NB,SB,P} +function BlockMultiFieldStyle(::BlockMultiFieldStyle{NB,SB,P},spaces) where {NB,SB,P} @check length(spaces) == sum(SB) return BlockMultiFieldStyle(NB,SB,P) end -function BlockMultiFieldStyle(::BlockMultiFieldStyle{0,0,0},spaces::Vector{<:SingleFieldFESpace}) +function BlockMultiFieldStyle(::BlockMultiFieldStyle{0,0,0},spaces) NB = length(spaces) return BlockMultiFieldStyle(NB) end diff --git a/src/ODEs/TransientFETools/TransientFESpaces.jl b/src/ODEs/TransientFETools/TransientFESpaces.jl index b13ebec12..b6312b96b 100644 --- a/src/ODEs/TransientFETools/TransientFESpaces.jl +++ b/src/ODEs/TransientFETools/TransientFESpaces.jl @@ -34,7 +34,7 @@ end Allocate the space to be used as first argument in evaluate! """ function allocate_trial_space(U::TransientTrialFESpace) - HomogeneousTrialFESpace(U.space) + HomogeneousTrialFESpace(U.space) end """ @@ -43,14 +43,14 @@ Time evaluation allocating Dirichlet vals function evaluate(U::TransientTrialFESpace,t::Real) Ut = allocate_trial_space(U) evaluate!(Ut,U,t) - Ut + return Ut end """ We can evaluate at `nothing` when we do not care about the Dirichlet vals """ function evaluate(U::TransientTrialFESpace,t::Nothing) - U.Ud0 + return U.Ud0 end evaluate(U::TrialFESpace,t::Nothing) = U @@ -80,6 +80,11 @@ Time 2nd derivative of the Dirichlet functions ∂tt(U::MultiFieldFESpace) = MultiFieldFESpace(∂tt.(U.spaces)) ∂tt(t::T) where T<:Number = zero(T) +zero_free_values(f::TransientTrialFESpace) = zero_free_values(f.space) +has_constraints(f::TransientTrialFESpace) = has_constraints(f.space) +get_dof_value_type(f::TransientTrialFESpace) = get_dof_value_type(f.space) +get_vector_type(f::TransientTrialFESpace) = get_vector_type(f.space) + # Testing the interface function test_transient_trial_fe_space(Uh) @@ -100,68 +105,122 @@ end # Define the TransientTrialFESpace interface for stationary spaces -function evaluate!(Ut::FESpace,U::FESpace,t::Real) - U -end +evaluate!(Ut::FESpace,U::FESpace,t::Real) = U +allocate_trial_space(U::FESpace) = U +evaluate(U::FESpace,t::Real) = U +evaluate(U::FESpace,t::Nothing) = U -function allocate_trial_space(U::FESpace) - U +@static if VERSION >= v"1.3" + (U::FESpace)(t) = U end -function evaluate(U::FESpace,t::Real) - U +# Define the interface for MultiField + +struct TransientMultiFieldTrialFESpace{MS<:MultiFieldStyle,CS<:ConstraintStyle,V} + vector_type::Type{V} + spaces::Vector + multi_field_style::MS + constraint_style::CS + function TransientMultiFieldTrialFESpace( + ::Type{V}, + spaces::Vector, + multi_field_style::MultiFieldStyle) where V + @assert length(spaces) > 0 + + MS = typeof(multi_field_style) + if any( map(has_constraints,spaces) ) + constraint_style = Constrained() + else + constraint_style = UnConstrained() + end + CS = typeof(constraint_style) + new{MS,CS,V}(V,spaces,multi_field_style,constraint_style) + end end -function evaluate(U::FESpace,t::Nothing) - U +# Default constructors +function TransientMultiFieldFESpace(spaces::Vector; + style = ConsecutiveMultiFieldStyle()) + Ts = map(get_dof_value_type,spaces) + T = typeof(*(map(zero,Ts)...)) + if isa(style,BlockMultiFieldStyle) + style = BlockMultiFieldStyle(style,spaces) + VT = typeof(mortar(map(zero_free_values,spaces))) + else + VT = Vector{T} + end + TransientMultiFieldTrialFESpace(VT,spaces,style) end -@static if VERSION >= v"1.3" - (U::FESpace)(t) = U +function TransientMultiFieldFESpace(::Type{V},spaces::Vector) where V + TransientMultiFieldTrialFESpace(V,spaces,ConsecutiveMultiFieldStyle()) end -# Define the interface for MultiField +function TransientMultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace}; + style = ConsecutiveMultiFieldStyle()) + MultiFieldFESpace(spaces,style=style) +end -struct TransientMultiFieldTrialFESpace - spaces::Vector +function TransientMultiFieldFESpace(::Type{V},spaces::Vector{<:SingleFieldFESpace}) where V + MultiFieldFESpace(V,spaces,ConsecutiveMultiFieldStyle()) end + Base.iterate(m::TransientMultiFieldTrialFESpace) = iterate(m.spaces) Base.iterate(m::TransientMultiFieldTrialFESpace,state) = iterate(m.spaces,state) Base.getindex(m::TransientMultiFieldTrialFESpace,field_id::Integer) = m.spaces[field_id] Base.length(m::TransientMultiFieldTrialFESpace) = length(m.spaces) - -function TransientMultiFieldFESpace(spaces::Vector) - TransientMultiFieldTrialFESpace(spaces) -end - -function TransientMultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace}) - MultiFieldFESpace(spaces) -end - function evaluate!(Ut::T,U::TransientMultiFieldTrialFESpace,t::Real) where T spaces_at_t = [evaluate!(Uti,Ui,t) for (Uti,Ui) in zip(Ut,U)] - MultiFieldFESpace(spaces_at_t) + mfs = MultiFieldStyle(U) + return MultiFieldFESpace(spaces_at_t;style=mfs) end function allocate_trial_space(U::TransientMultiFieldTrialFESpace) spaces = allocate_trial_space.(U.spaces) - MultiFieldFESpace(spaces) + mfs = MultiFieldStyle(U) + return MultiFieldFESpace(spaces;style=mfs) end function evaluate(U::TransientMultiFieldTrialFESpace,t::Real) Ut = allocate_trial_space(U) evaluate!(Ut,U,t) - Ut + return Ut end function evaluate(U::TransientMultiFieldTrialFESpace,t::Nothing) - MultiFieldFESpace([evaluate(fesp,nothing) for fesp in U.spaces]) + spaces = [evaluate(fesp,nothing) for fesp in U.spaces] + mfs = MultiFieldStyle(U) + MultiFieldFESpace(spaces;style=mfs) end (U::TransientMultiFieldTrialFESpace)(t) = evaluate(U,t) function ∂t(U::TransientMultiFieldTrialFESpace) spaces = ∂t.(U.spaces) - TransientMultiFieldFESpace(spaces) + mfs = MultiFieldStyle(U) + TransientMultiFieldFESpace(spaces;style=mfs) +end + +function zero_free_values(f::TransientMultiFieldTrialFESpace{<:BlockMultiFieldStyle{NB,SB,P}}) where {NB,SB,P} + block_ranges = get_block_ranges(NB,SB,P) + block_num_dofs = map(range->sum(map(num_free_dofs,f.spaces[range])),block_ranges) + block_vtypes = map(range->get_vector_type(first(f.spaces[range])),block_ranges) + return mortar(map(allocate_vector,block_vtypes,block_num_dofs)) +end + +get_dof_value_type(f::TransientMultiFieldTrialFESpace{MS,CS,V}) where {MS,CS,V} = eltype(V) +get_vector_type(f::TransientMultiFieldTrialFESpace) = f.vector_type +ConstraintStyle(::Type{TransientMultiFieldTrialFESpace{S,B,V}}) where {S,B,V} = B() +ConstraintStyle(::TransientMultiFieldTrialFESpace) = ConstraintStyle(typeof(f)) +MultiFieldStyle(::Type{TransientMultiFieldTrialFESpace{S,B,V}}) where {S,B,V} = S() +MultiFieldStyle(f::TransientMultiFieldTrialFESpace) = MultiFieldStyle(typeof(f)) + +function SparseMatrixAssembler(mat,vec, + trial::TransientMultiFieldTrialFESpace{MS}, + test ::TransientMultiFieldTrialFESpace{MS}, + strategy::AssemblyStrategy=DefaultAssemblyStrategy() + ) where MS <: BlockMultiFieldStyle + mfs = MultiFieldStyle(test) + return BlockSparseMatrixAssembler(mfs,trial,test,SparseMatrixBuilder(mat),ArrayBuilder(vec),strategy) end diff --git a/src/ODEs/TransientFETools/TransientFETools.jl b/src/ODEs/TransientFETools/TransientFETools.jl index 2b262b477..08f862407 100644 --- a/src/ODEs/TransientFETools/TransientFETools.jl +++ b/src/ODEs/TransientFETools/TransientFETools.jl @@ -112,6 +112,13 @@ import Gridap.CellData: gradient import Gridap.CellData: ∇∇ import Gridap.CellData: change_domain import Gridap.FESpaces: BasisStyle +using Gridap.FESpaces: Constrained, UnConstrained, AssemblyStrategy +using Gridap.MultiField: ConsecutiveMultiFieldStyle, BlockSparseMatrixAssembler +import Gridap.MultiField: ConstraintStyle, MultiFieldStyle, BlockMultiFieldStyle +import Gridap.FESpaces: zero_free_values, has_constraints, SparseMatrixAssembler +import Gridap.FESpaces: get_dof_value_type, get_vector_type + +using BlockArrays include("TransientFESpaces.jl") diff --git a/test/ODEsTests/TransientFEsTests/TransientBlockMultiFieldStyleTests.jl b/test/ODEsTests/TransientFEsTests/TransientBlockMultiFieldStyleTests.jl new file mode 100644 index 000000000..6c6e4dbe3 --- /dev/null +++ b/test/ODEsTests/TransientFEsTests/TransientBlockMultiFieldStyleTests.jl @@ -0,0 +1,103 @@ +module TransientBlockMultiFieldStyleTests +using Test, BlockArrays, SparseArrays, LinearAlgebra + +using Gridap +using Gridap.FESpaces, Gridap.ReferenceFEs, Gridap.MultiField +using Gridap.ODEs.TransientFETools + +function main(n_spaces,mfs,weakform,Ω,dΩ,U,V) + mass, biform, liform = weakform + res(t,x,y) = mass(t,∂t(x),y) + biform(t,x,y) - liform(t,y) + jac(t,x,dx,y) = biform(t,dx,y) + jac_t(t,xt,dxt,y) = mass(t,dxt,y) + + ############################################################################################ + # Normal assembly + + Y = MultiFieldFESpace(fill(V,n_spaces)) + X = TransientMultiFieldFESpace(fill(U,n_spaces)) + + u = get_trial_fe_basis(X(0.0)) + v = get_fe_basis(Y) + uₜ = TransientCellField(u,(u,)) + + matdata_jac = collect_cell_matrix(X(0),Y,jac(0,uₜ,u,v)) + matdata_jac_t = collect_cell_matrix(X(0),Y,jac_t(0,uₜ,u,v)) + matdata_jacs = (matdata_jac,matdata_jac_t) + matdata = TransientFETools._vcat_matdata(matdata_jacs) + vecdata = collect_cell_vector(Y,liform(0,v)) + + assem = SparseMatrixAssembler(X(0),Y) + A1 = assemble_matrix(assem,matdata) + b1 = assemble_vector(assem,vecdata) + + ############################################################################################ + # Block MultiFieldStyle + + Yb = MultiFieldFESpace(fill(V,n_spaces);style=mfs) + Xb = TransientMultiFieldFESpace(fill(U,n_spaces);style=mfs) + test_fe_space(Yb) + test_fe_space(Xb(0)) + + ub = get_trial_fe_basis(Xb(0)) + vb = get_fe_basis(Yb) + ubₜ = TransientCellField(ub,(ub,)) + + bmatdata_jac = collect_cell_matrix(Xb(0),Yb,jac(0,ubₜ,ub,vb)) + bmatdata_jac_t = collect_cell_matrix(Xb(0),Yb,jac_t(0,ubₜ,ub,vb)) + bmatdata_jacs = (bmatdata_jac,bmatdata_jac_t) + bmatdata = TransientFETools._vcat_matdata(bmatdata_jacs) + bvecdata = collect_cell_vector(Yb,liform(0,vb)) + + ############################################################################################ + # Block Assembly + + assem_blocks = SparseMatrixAssembler(Xb,Yb) + + A1_blocks = assemble_matrix(assem_blocks,bmatdata) + b1_blocks = assemble_vector(assem_blocks,bvecdata) + @test A1 ≈ A1_blocks + @test b1 ≈ b1_blocks + + y1_blocks = similar(b1_blocks) + mul!(y1_blocks,A1_blocks,b1_blocks) + y1 = similar(b1) + mul!(y1,A1,b1) + @test y1_blocks ≈ y1 + + A3_blocks = allocate_matrix(assem_blocks,bmatdata) + b3_blocks = allocate_vector(assem_blocks,bvecdata) + assemble_matrix!(A3_blocks,assem_blocks,bmatdata) + assemble_vector!(b3_blocks,assem_blocks,bvecdata) + @test A3_blocks ≈ A1 + @test b3_blocks ≈ b1_blocks + +end + +############################################################################################ + +sol(x,t) = sum(x) +sol(t::Real) = x->sol(x,t) + +model = CartesianDiscreteModel((0.0,1.0,0.0,1.0),(5,5)) +Ω = Triangulation(model) + +reffe = LagrangianRefFE(Float64,QUAD,1) +V = FESpace(Ω, reffe; dirichlet_tags="boundary") +U = TransientTrialFESpace(V,sol) + +dΩ = Measure(Ω, 2) +mass2(t,(u1t,u2t),(v1,v2)) = ∫(u1t⋅v1)*dΩ +biform2(t,(u1,u2),(v1,v2)) = ∫(∇(u1)⋅∇(v1) + u2⋅v2 - u1⋅v2)*dΩ +liform2(t,(v1,v2)) = ∫(v1 - v2)*dΩ +mass3(t,(u1t,u2t,u3t),(v1,v2,v3)) = ∫(u1t⋅v1)*dΩ +biform3(t,(u1,u2,u3),(v1,v2,v3)) = ∫(∇(u1)⋅∇(v1) + u2⋅v2 - u1⋅v2 - u3⋅v2 - u2⋅v3)*dΩ +liform3(t,(v1,v2,v3)) = ∫(v1 - v2 + 2.0*v3)*dΩ + +for (n_spaces,weakform) in zip([2,3],[(mass2,biform2,liform2),(mass3,biform3,liform3)]) + for mfs in [BlockMultiFieldStyle(),BlockMultiFieldStyle(2,(1,n_spaces-1))] + main(n_spaces,mfs,weakform,Ω,dΩ,U,V) + end +end + +end # module