Merge pull request #946 from gridap/Transient_MultifieldStyle

Transient multifield style

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

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

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

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

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")
 

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