Merge pull request #185 from gridap/create_trial_fe_space_out_of_distributed_cell_field

Necessary changes required to be able to instantiate a TrialFESpace with an instance of DistributedCellField

Author: amartinhuertas

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 
+
+### Added
+
+- New overloads for the `TrialFESpace` constructor where the data to be imposed is passed as a `DistributedCellField`. Since PR[#183](https://github.com/gridap/GridapDistributed.jl/pull/183).
+
 ## [0.4.10] - 2025-09-29
 
 ### Added

Project.toml

@@ -21,7 +21,7 @@ BlockArrays = "1"
 CircularArrays = "1.4.0"
 FillArrays = "1"
 ForwardDiff = "0.10, 1"
-Gridap = "0.19.4"
+Gridap = "0.19.6"
 LinearAlgebra = "1"
 MPI = "0.16, 0.17, 0.18, 0.19, 0.20"
 PartitionedArrays = "0.3.3"

src/FESpaces.jl

@@ -414,6 +414,21 @@ function FESpaces.TrialFESpace!(f::DistributedSingleFieldFESpace,fun)
   DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type,f.metadata)
 end
 
+function FESpaces.TrialFESpace(f::DistributedSingleFieldFESpace,cf::DistributedCellField)
+  spaces = map(local_views(f),local_views(cf)) do s, field
+    TrialFESpace(s,field)
+  end
+  DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type,f.metadata)
+end
+
+function FESpaces.TrialFESpace(cf::DistributedCellField,f::DistributedSingleFieldFESpace)
+  spaces = map(local_views(f),local_views(cf)) do s, field
+    TrialFESpace(s,field)
+  end
+  DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type,f.metadata)
+end
+
+
 function FESpaces.HomogeneousTrialFESpace(f::DistributedSingleFieldFESpace)
   spaces = map(f.spaces) do s
     HomogeneousTrialFESpace(s)
@@ -478,6 +493,16 @@ function FESpaces.interpolate_dirichlet!(
   FEFunction(f,free_values,dirichlet_values)
 end
 
+function FESpaces.interpolate_dirichlet!(
+  u::DistributedCellField, free_values::AbstractVector,
+  dirichlet_values::AbstractArray{<:AbstractVector},
+  f::DistributedSingleFieldFESpace)
+  map(local_views(u), f.spaces,local_views(free_values),dirichlet_values) do u,V,fvec,dvec
+    interpolate_dirichlet!(u,fvec,dvec,V)
+  end
+  FEFunction(f,free_values,dirichlet_values)
+end
+
 function FESpaces.interpolate_everywhere(u, f::DistributedSingleFieldFESpace)
   free_values = zero_free_values(f)
   dirichlet_values = get_dirichlet_dof_values(f)

test/FESpacesTests.jl

@@ -85,66 +85,87 @@ function main(distribute,parts,das)
   model = CartesianDiscreteModel(ranks,parts,domain,cells)
   Ω = Triangulation(model)
   Γ = Boundary(model)
+  dΩ = Measure(Ω,3)
+
+  function _interpolation_tests(u,model,Ω,Γ,output)
+      reffe = ReferenceFE(lagrangian,Float64,1)
+      V = TestFESpace(model,reffe,dirichlet_tags="boundary")
+      U = TrialFESpace(u,V)
+      V2 = FESpace(Ω,reffe)
+      @test get_vector_type(V) <: PVector{<:Vector}
+      @test get_vector_type(U) <: PVector{<:Vector}
+      @test get_vector_type(V2) <: PVector{<:Vector}
+
+      free_values_partition = map(partition(V.gids)) do indices 
+        ones(Float64,local_length(indices))
+      end 
+
+      free_values = PVector(free_values_partition,partition(V.gids))
+      fh = FEFunction(U,free_values)
+      zh = zero(U)
+      uh = interpolate(u,U)
+      eh = u - uh
+
+      uh_dir = interpolate_dirichlet(u,U)
+      free_values = zero_free_values(U)
+      dirichlet_values = get_dirichlet_dof_values(U)
+      uh_dir2 = interpolate_dirichlet!(u,free_values,dirichlet_values,U)
+
+      uh_everywhere = interpolate_everywhere(u,U)
+      dirichlet_values0 = zero_dirichlet_values(U)
+      uh_everywhere_ = interpolate_everywhere!(u,free_values,dirichlet_values0,U)
+      eh2 = u - uh_everywhere
+      eh2_ = u - uh_everywhere_
+
+      uh_everywhere2 = interpolate_everywhere(uh_everywhere,U)
+      uh_everywhere2_ = interpolate_everywhere!(uh_everywhere,free_values,dirichlet_values,U)
+      eh3 = u - uh_everywhere2
+
+      dofs      = get_fe_dof_basis(U)
+      cell_vals = dofs(uh)
+      gather_free_values!(free_values,U,cell_vals)
+      gather_free_and_dirichlet_values!(free_values,dirichlet_values,U,cell_vals)
+      uh4 = FEFunction(U,free_values,dirichlet_values)
+      eh4 = u - uh4
+
+      dΩ = Measure(Ω,3)
+      cont   = ∫( abs2(eh) )dΩ
+      cont2  = ∫( abs2(eh2) )dΩ
+      cont2_ = ∫( abs2(eh2_) )dΩ
+      cont3  = ∫( abs2(eh3) )dΩ
+      cont4  = ∫( abs2(eh4) )dΩ
+      @test sqrt(sum(cont))   < 1.0e-9
+      @test sqrt(sum(cont2))  < 1.0e-9
+      @test sqrt(sum(cont2_)) < 1.0e-9
+      @test sqrt(sum(cont3))  < 1.0e-9
+      @test sqrt(sum(cont4))  < 1.0e-9
+
+      writevtk(Ω,joinpath(output,"Ω"), nsubcells=10,
+               celldata=["err"=>cont[Ω]],
+               cellfields=["uh"=>uh,"zh"=>zh,"eh"=>eh])
+
+      writevtk(Γ,joinpath(output,"Γ"),cellfields=["uh"=>uh])
+  end
 
+  # Testing interpolation with a Julia function
   u((x,y)) = x+y
+  _interpolation_tests(u,model,Ω,Γ,output)
+
+  # Testing interpolation with a DistributedCellField
+  # The DistributedCellField MUST be defined in all cells
+  # (i.e., owned + ghost) for the tests within 
+  # _interpolation_tests(...) to pass. The problematic part is 
+  # in the interpolation of the dirichlet values
+  Ωghosts = Triangulation(with_ghost, model)
+  fields = map(local_views(Ωghosts)) do Ω
+    CellField(u,Ω)
+  end
+  u_cf = GridapDistributed.DistributedCellField(fields,Ωghosts)
+  _interpolation_tests(u_cf,model,Ωghosts,Γ,output)
+
   reffe = ReferenceFE(lagrangian,Float64,1)
   V = TestFESpace(model,reffe,dirichlet_tags="boundary")
   U = TrialFESpace(u,V)
-  V2 = FESpace(Ω,reffe)
-  @test get_vector_type(V) <: PVector{<:Vector}
-  @test get_vector_type(U) <: PVector{<:Vector}
-  @test get_vector_type(V2) <: PVector{<:Vector}
-
-  free_values_partition = map(partition(V.gids)) do indices 
-    ones(Float64,local_length(indices))
-  end 
-
-  free_values = PVector(free_values_partition,partition(V.gids))
-  fh = FEFunction(U,free_values)
-  zh = zero(U)
-  uh = interpolate(u,U)
-  eh = u - uh
-
-  uh_dir = interpolate_dirichlet(u,U)
-  free_values = zero_free_values(U)
-  dirichlet_values = get_dirichlet_dof_values(U)
-  uh_dir2 = interpolate_dirichlet!(u,free_values,dirichlet_values,U)
-
-  uh_everywhere = interpolate_everywhere(u,U)
-  dirichlet_values0 = zero_dirichlet_values(U)
-  uh_everywhere_ = interpolate_everywhere!(u,free_values,dirichlet_values0,U)
-  eh2 = u - uh_everywhere
-  eh2_ = u - uh_everywhere_
-
-  uh_everywhere2 = interpolate_everywhere(uh_everywhere,U)
-  uh_everywhere2_ = interpolate_everywhere!(uh_everywhere,free_values,dirichlet_values,U)
-  eh3 = u - uh_everywhere2
-
-  dofs      = get_fe_dof_basis(U)
-  cell_vals = dofs(uh)
-  gather_free_values!(free_values,U,cell_vals)
-  gather_free_and_dirichlet_values!(free_values,dirichlet_values,U,cell_vals)
-  uh4 = FEFunction(U,free_values,dirichlet_values)
-  eh4 = u - uh4
-
-  dΩ = Measure(Ω,3)
-  cont   = ∫( abs2(eh) )dΩ
-  cont2  = ∫( abs2(eh2) )dΩ
-  cont2_ = ∫( abs2(eh2_) )dΩ
-  cont3  = ∫( abs2(eh3) )dΩ
-  cont4  = ∫( abs2(eh4) )dΩ
-  @test sqrt(sum(cont))   < 1.0e-9
-  @test sqrt(sum(cont2))  < 1.0e-9
-  @test sqrt(sum(cont2_)) < 1.0e-9
-  @test sqrt(sum(cont3))  < 1.0e-9
-  @test sqrt(sum(cont4))  < 1.0e-9
-
-
-  writevtk(Ω,joinpath(output,"Ω"), nsubcells=10,
-           celldata=["err"=>cont[Ω]],
-           cellfields=["uh"=>uh,"zh"=>zh,"eh"=>eh])
-
-  writevtk(Γ,joinpath(output,"Γ"),cellfields=["uh"=>uh])
 
   # Assembly
   Ωass  = Triangulation(das,model)