Interpolate multifield

NEWS.md

@@ -12,8 +12,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   - Operator `⊗` (\otimes) as an alias of `outer`. Since PR [#239](/~https://github.com/gridap/Gridap.jl/pull/239).
   - Support for (symmetric) 4th order tensors. Since PR [#239](/~https://github.com/gridap/Gridap.jl/pull/239).
   - Optimizations for symmetric 2nd order tensors. Since PR [#239](/~https://github.com/gridap/Gridap.jl/pull/239).
-  - Methods for `cross` function (aka `×` (\times)) to operate with `VectorValues`.
-  Since PR [#280](/~https://github.com/gridap/Gridap.jl/pull/280).
+  - Methods for `cross` function (aka `×` (\times)) to operate with `VectorValues`. Since PR [#280](/~https://github.com/gridap/Gridap.jl/pull/280).
+  - Interpolation is now supported also for multifield spaces. Since PR [#279](/~https://github.com/gridap/Gridap.jl/pull/279).
 
 ### Changed
 
@@ -26,6 +26,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
     - Gridap re-exports `dot`, `⋅`, and other names from LinearAlbegra that are used 
   often in Gridap code. 
     - Function `n_components` is renamed to `num_components`.
+  - The `SingleFieldFESpace` interface has changed. The function `gather_free_and_dirichlet_values!` 
+  has been added as mandatory for all FE space implementations and the old function `gather_free_and_dirichlet_values` 
+  is now optional. Since PR [#279](/~https://github.com/gridap/Gridap.jl/pull/279).
 
 ## [0.10.4] - 2020-6-8
 

src/FESpaces/CLagrangianFESpaces.jl

@@ -71,12 +71,12 @@ function scatter_free_and_dirichlet_values(f::CLagrangianFESpace,free_values,dir
   scatter_free_and_dirichlet_values(f.space,free_values,dirichlet_values)
 end
 
-function gather_free_and_dirichlet_values(f::CLagrangianFESpace,cell_vals)
-  gather_free_and_dirichlet_values(f.space,cell_vals)
+function gather_free_and_dirichlet_values!(free_values,dirichlet_values,f::CLagrangianFESpace,cell_vals)
+  gather_free_and_dirichlet_values!(free_values,dirichlet_values,f.space,cell_vals)
 end
 
-function gather_dirichlet_values(f::CLagrangianFESpace,cell_vals)
-  gather_dirichlet_values(f.space,cell_vals)
+function gather_dirichlet_values!(dirichlet_vals,f::CLagrangianFESpace,cell_vals)
+    gather_dirichlet_values!(dirichlet_vals,f.space,cell_vals)
 end
 
 # Helpers

src/FESpaces/DirichletFESpaces.jl

@@ -60,9 +60,9 @@ function scatter_free_and_dirichlet_values(f::DirichletFESpace,fv,dv)
   scatter_free_and_dirichlet_values(f.space,dv,fv)
 end
 
-function gather_free_and_dirichlet_values(f::DirichletFESpace,cv)
-  dv, fv = gather_free_and_dirichlet_values(f.space,cv)
-  (fv, dv)
+function gather_free_and_dirichlet_values!(fv,dv,f::DirichletFESpace,cv)
+  gather_free_and_dirichlet_values!(dv,fv,f.space,cv)
+  (fv,dv)
 end
 
 function TrialFESpace(f::DirichletFESpace)

src/FESpaces/ExtendedFESpaces.jl

@@ -225,9 +225,9 @@ function scatter_free_and_dirichlet_values(f::ExtendedFESpace,fv,dv)
     f.trian.cell_to_oldcell)
 end
 
-function gather_free_and_dirichlet_values(f::ExtendedFESpace,cv)
+function gather_free_and_dirichlet_values!(fv,dv,f::ExtendedFESpace,cv)
   _cv = reindex(cv,f.trian.cell_to_oldcell)
-  gather_free_and_dirichlet_values(f.space,_cv)
+  gather_free_and_dirichlet_values!(fv,dv,f.space,_cv)
 end
 
 # Delegated functions

src/FESpaces/FESpaces.jl

@@ -141,11 +141,14 @@ export num_dirichlet_dofs
 export get_cell_dofs
 export zero_dirichlet_values
 export gather_free_and_dirichlet_values
+export gather_free_and_dirichlet_values!
 export scatter_free_and_dirichlet_values
 export get_dirichlet_values
 export gather_dirichlet_values
+export gather_dirichlet_values!
 export num_dirichlet_tags
 export gather_free_values
+export gather_free_values!
 export get_dirichlet_dof_tag
 export compute_free_and_dirichlet_values
 export compute_dirichlet_values

src/FESpaces/FESpacesWithLastDofRemoved.jl

@@ -64,6 +64,16 @@ function gather_free_and_dirichlet_values(
   (fv, dv)
 end
 
+function gather_free_and_dirichlet_values!(
+  fv,dv,f::FESpaceWithLastDofRemoved,cv)
+  _fv , _dv = gather_free_and_dirichlet_values(f.space,cv)
+  @assert length(_dv) == 0
+  l = length(_fv)
+  fv .= SubVector(_fv,1,l-1)
+  dv .= SubVector(_fv,l,l)
+  (fv, dv)
+end
+
 function TrialFESpace(f::FESpaceWithLastDofRemoved)
   U = TrialFESpace(f.space)
   FESpaceWithLastDofRemoved(U)

src/FESpaces/FESpacesWithLinearConstraints.jl

@@ -419,6 +419,19 @@ function gather_free_and_dirichlet_values(f::FESpaceWithLinearConstraints,cell_t
   fmdof_to_val, dmdof_to_val
 end
 
+function gather_free_and_dirichlet_values!(fmdof_to_val,dmdof_to_val,f::FESpaceWithLinearConstraints,cell_to_ludof_to_val)
+  fdof_to_val, ddof_to_val = gather_free_and_dirichlet_values(f.space,cell_to_ludof_to_val)
+  _setup_mdof_to_val!(
+    fmdof_to_val,
+    dmdof_to_val,
+    fdof_to_val,
+    ddof_to_val,
+    f.mDOF_to_DOF,
+    f.n_fdofs,
+    f.n_fmdofs)
+  fmdof_to_val, dmdof_to_val
+end
+
 function _setup_mdof_to_val!(
   fmdof_to_val,
   dmdof_to_val,

src/FESpaces/SingleFieldFESpaces.jl

@@ -47,11 +47,10 @@ end
 
 """
 """
-function gather_free_and_dirichlet_values(f::SingleFieldFESpace,cell_vals)
+function gather_free_and_dirichlet_values!(free_values, dirichlet_values,fs::SingleFieldFESpace,cell_vals)
   @abstractmethod
 end
 
-
 """
 """
 function test_single_field_fe_space(f::SingleFieldFESpace,pred=(==))
@@ -66,6 +65,12 @@ function test_single_field_fe_space(f::SingleFieldFESpace,pred=(==))
   fv, dv = gather_free_and_dirichlet_values(f,cell_vals)
   @test pred(fv,free_values)
   @test pred(dv,dirichlet_values)
+  gather_free_and_dirichlet_values!(fv,dv,f,cell_vals)
+  @test pred(fv,free_values)
+  @test pred(dv,dirichlet_values)
+  fv, dv = gather_free_and_dirichlet_values!(fv,dv,f,cell_vals)
+  @test pred(fv,free_values)
+  @test pred(dv,dirichlet_values)
   fe_function = FEFunction(f,free_values,dirichlet_values)
   @test isa(fe_function, SingleFieldFEFunction)
   test_fe_function(fe_function)
@@ -237,28 +242,61 @@ function get_dirichlet_values(f::SingleFieldFESpace)
   zero_dirichlet_values(f)
 end
 
+"""
+"""
+function gather_free_and_dirichlet_values(f::SingleFieldFESpace,cell_vals)
+  free_values = zero_free_values(f)
+  dirichlet_values = zero_dirichlet_values(f)
+  gather_free_and_dirichlet_values!(free_values,dirichlet_values,f,cell_vals)
+end
+
 """
 """
 function gather_dirichlet_values(f::SingleFieldFESpace,cell_vals)
-  _, dirichlet_values = gather_free_and_dirichlet_values(f,cell_vals)
+  dirichlet_values = zero_dirichlet_values(f)
+  gather_dirichlet_values!(dirichlet_values,f,cell_vals)
   dirichlet_values
 end
 
 """
 """
 function gather_free_values(f::SingleFieldFESpace,cell_vals)
-  free_values, _ = gather_free_and_dirichlet_values(f,cell_vals)
+  free_values = zero_free_values(f)
+  gather_free_values!(free_values,f,cell_vals)
   free_values
 end
 
+"""
+"""
+function gather_dirichlet_values!(dirichlet_values,f::SingleFieldFESpace,cell_vals)
+  free_values = zero_free_values(f)
+  gather_free_and_dirichlet_values!(free_values,dirichlet_values,f,cell_vals)
+  dirichlet_values
+end
+
+"""
+"""
+function gather_free_values!(free_values,f::SingleFieldFESpace,cell_vals)
+    dirichlet_values = zero_dirichlet_values(f)
+    gather_free_and_dirichlet_values!(free_values,dirichlet_values,f,cell_vals)
+    free_values
+end
+
 """
 The resulting FE function is in the space (in particular it fulfills Dirichlet BCs
 even in the case that the given cell field does not fulfill them)
 """
 function interpolate(fs::SingleFieldFESpace,object)
-  cell_vals = _cell_vals(fs,object)
-  free_values = gather_free_values(fs,cell_vals)
-  FEFunction(fs,free_values)
+  free_values = zero_free_values(fs)
+  interpolate!(free_values,fs,object)
+end
+
+"""
+"""
+function interpolate!(free_values,fs::SingleFieldFESpace,object)
+    cell_vals = _cell_vals(fs,object)
+    gather_free_values!(free_values,fs,cell_vals)
+    FEFunction(fs,free_values)
 end
 
 function _cell_vals(fs::SingleFieldFESpace,object)
@@ -273,17 +311,33 @@ like interpolate, but also compute new degrees of freedom for the dirichlet comp
 The resulting FEFunction does not necessary belongs to the underlying space
 """
 function interpolate_everywhere(fs::SingleFieldFESpace,object)
+  free_values = zero_free_values(fs)
+  dirichlet_values = zero_dirichlet_values(fs)
+  interpolate_everywhere!(free_values,dirichlet_values,fs,object)
+end
+
+"""
+"""
+function interpolate_everywhere!(free_values,dirichlet_values,fs::SingleFieldFESpace,object)
   cell_vals = _cell_vals(fs,object)
-  free_values, dirichlet_values = gather_free_and_dirichlet_values(fs,cell_vals)
+  gather_free_and_dirichlet_values!(free_values,dirichlet_values,fs,cell_vals)
   FEFunction(fs,free_values,dirichlet_values)
 end
 
 """
 """
 function interpolate_dirichlet(fs::SingleFieldFESpace,object)
-  cell_vals = _cell_vals(fs,object)
-  dirichlet_values = gather_dirichlet_values(fs,cell_vals)
   free_values = zero_free_values(fs)
+  dirichlet_values = zero_dirichlet_values(fs)
+  interpolate_dirichlet!(free_values,dirichlet_values,fs,object)
+end
+
+"""
+"""
+function interpolate_dirichlet!(free_values,dirichlet_values,fs::SingleFieldFESpace,object)
+  cell_vals = _cell_vals(fs,object)
+  gather_dirichlet_values!(dirichlet_values,fs,cell_vals)
+  fill!(free_values,zero(eltype(free_values)))
   FEFunction(fs,free_values, dirichlet_values)
 end
 

src/FESpaces/TrialFESpaces.jl

@@ -94,10 +94,16 @@ scatter_free_and_dirichlet_values(f::TrialFESpace,fv,dv) = scatter_free_and_diri
 
 gather_free_and_dirichlet_values(f::TrialFESpace,cv) = gather_free_and_dirichlet_values(f.space,cv)
 
+gather_free_and_dirichlet_values!(fv,dv,f::TrialFESpace,cv) = gather_free_and_dirichlet_values!(fv,dv,f.space,cv)
+
 gather_dirichlet_values(f::TrialFESpace,cv) = gather_dirichlet_values(f.space,cv)
 
+gather_dirichlet_values!(dv,f::TrialFESpace,cv) = gather_dirichlet_values!(dv,f.space,cv)
+
 gather_free_values(f::TrialFESpace,cv) = gather_free_values(f.space,cv)
 
+gather_free_values!(fv,f::TrialFESpace,cv) = gather_free_values!(fv,f.space,cv)
+
 function get_constraint_kernel_matrix_cols(f::TrialFESpace)
   get_constraint_kernel_matrix_cols(f.space)
 end

src/FESpaces/UnconstrainedFESpaces.jl

@@ -89,13 +89,11 @@ function scatter_free_and_dirichlet_values(f::UnconstrainedFESpace,free_values,d
   LocalToGlobalPosNegArray(cell_dofs,free_values,dirichlet_values)
 end
 
-function gather_free_and_dirichlet_values(f::UnconstrainedFESpace,cell_vals)
+function gather_free_and_dirichlet_values!(free_vals,dirichlet_vals,f::UnconstrainedFESpace,cell_vals)
 
   cell_dofs = get_cell_dofs(f)
   cache_vals = array_cache(cell_vals)
   cache_dofs = array_cache(cell_dofs)
-  free_vals = zero_free_values(f)
-  dirichlet_vals = zero_dirichlet_values(f)
   cells = 1:length(cell_vals)
 
   _free_and_dirichlet_values_fill!(
@@ -107,16 +105,15 @@ function gather_free_and_dirichlet_values(f::UnconstrainedFESpace,cell_vals)
     cell_dofs,
     cells)
 
-  (free_vals, dirichlet_vals)
+  (free_vals,dirichlet_vals)
 end
 
-function gather_dirichlet_values(f::UnconstrainedFESpace,cell_vals)
+function gather_dirichlet_values!(dirichlet_vals,f::UnconstrainedFESpace,cell_vals)
 
   cell_dofs = get_cell_dofs(f)
   cache_vals = array_cache(cell_vals)
   cache_dofs = array_cache(cell_dofs)
   free_vals = zero_free_values(f)
-  dirichlet_vals = zero_dirichlet_values(f)
   cells = f.dirichlet_cells
 
   _free_and_dirichlet_values_fill!(

src/FESpaces/ZeroMeanFESpaces.jl

@@ -112,7 +112,13 @@ scatter_free_and_dirichlet_values(f::ZeroMeanFESpace,fv,dv) = scatter_free_and_d
 
 gather_free_and_dirichlet_values(f::ZeroMeanFESpace,cv) = gather_free_and_dirichlet_values(f.space,cv)
 
+gather_free_and_dirichlet_values!(fv,dv,f::ZeroMeanFESpace,cv) = gather_free_and_dirichlet_values!(fv,dv,f.space,cv)
+
 gather_dirichlet_values(f::ZeroMeanFESpace,cv) = gather_dirichlet_values(f.space,cv)
 
+gather_dirichlet_values!(dv,f::ZeroMeanFESpace,cv) = gather_dirichlet_values!(dv,f.space,cv)
+
 gather_free_values(f::ZeroMeanFESpace,cv) = gather_free_values(f.space,cv)
 
+gather_free_values!(fv,f::ZeroMeanFESpace,cv) = gather_free_values!(fv,f.space,cv)
+

src/MultiField/MultiField.jl

@@ -23,6 +23,9 @@ using Gridap.Arrays: IdentityVector
 using Gridap.FESpaces: _operate_cell_basis
 using Gridap.FESpaces: _operate_cell_matrix_field
 using Gridap.FESpaces: SkeletonCellBasis
+using Gridap.FESpaces: interpolate!
+using Gridap.FESpaces: interpolate_everywhere!
+using Gridap.FESpaces: interpolate_dirichlet!
 
 import Gridap.Helpers: operate
 import Gridap.Arrays: get_array
@@ -86,6 +89,9 @@ import Gridap.FESpaces: fill_matrix_coo_numeric!
 import Gridap.FESpaces: fill_matrix_and_vector_coo_numeric!
 import Gridap.FESpaces: count_matrix_and_vector_nnz_coo
 import Gridap.FESpaces: fill_matrix_and_vector_coo_symbolic!
+import Gridap.FESpaces: interpolate
+import Gridap.FESpaces: interpolate_everywhere
+import Gridap.FESpaces: interpolate_dirichlet
 
 import Base: +, -
 

src/MultiField/MultiFieldFESpaces.jl

@@ -309,3 +309,48 @@ function compute_field_offsets(f::MultiFieldFESpace)
   offsets
 end
 
+"""
+The resulting MultiFieldFEFunction is in the space (in particular it fulfills Dirichlet BCs
+even in the case that the given cell field does not fulfill them)
+"""
+function interpolate(fe::MultiFieldFESpace,objects)
+  free_values = zero_free_values(fe)
+  blocks = SingleFieldFEFunction[]
+  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
+    free_values_i = restrict_to_field(fe,free_values,field)
+    uhi = interpolate!(free_values_i,U,object)
+    push!(blocks,uhi)
+  end
+  MultiFieldFEFunction(free_values,fe,blocks)
+end
+
+"""
+like interpolate, but also compute new degrees of freedom for the dirichlet component.
+The resulting MultiFieldFEFunction does not necessary belongs to the underlying space
+"""
+function interpolate_everywhere(fe::MultiFieldFESpace,objects)
+  free_values = zero_free_values(fe)
+  blocks = SingleFieldFEFunction[]
+  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
+    free_values_i = restrict_to_field(fe,free_values,field)
+    dirichlet_values_i = zero_dirichlet_values(U)
+    uhi = interpolate_everywhere!(free_values_i,dirichlet_values_i,U,object)
+    push!(blocks,uhi)
+  end
+  MultiFieldFEFunction(free_values,fe,blocks)
+end
+
+"""
+"""
+function interpolate_dirichlet(fe::MultiFieldFESpace,objects)
+  free_values = zero_free_values(fe)
+  blocks = SingleFieldFEFunction[]
+  for (field, (U,object)) in enumerate(zip(fe.spaces,objects))
+    free_values_i = restrict_to_field(fe,free_values,field)
+    dirichlet_values_i = zero_dirichlet_values(U)
+    uhi = interpolate_dirichlet!(free_values_i,dirichlet_values_i,U,object)
+    push!(blocks,uhi)
+  end
+  MultiFieldFEFunction(free_values,fe,blocks)
+end
+

test/MultiFieldTests/MultiFieldFESpacesTests.jl

@@ -85,6 +85,11 @@ for ref_st in ref_style
   cell_dofs_new = reindex(cell_dofs,cellids)
   @test isa(cell_dofs_new,MultiFieldCellArray)
   @test cell_dofs_new.block_ids == cell_dofs.block_ids
+
+  f(x) = sin(4*pi*(x[1]-x[2]^2))+1
+  fh = interpolate(X,[f,f])
+  fh = interpolate_everywhere(X,[f,f])
+  fh = interpolate_dirichlet(X,[f,f])
 end
 
 end # module