From 0e6cb126f9df9287313344bc1487085479f692af Mon Sep 17 00:00:00 2001 From: "Alberto F. Martin" Date: Sun, 2 Nov 2025 17:16:44 +1100 Subject: [PATCH 1/2] Necessary changes required to be able to instantiate a TrialFESpace with an instance of DistributedCellField --- NEWS.md | 6 ++ src/FESpaces.jl | 25 ++++++++ test/FESpacesTests.jl | 131 ++++++++++++++++++++++++------------------ 3 files changed, 107 insertions(+), 55 deletions(-) diff --git a/NEWS.md b/NEWS.md index c2e2b11d..b1eb72e9 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 + +### Added + +- New overloads for the `TrialFESpace` constructor where the data to be imposed is passed as a `DistributedCellField`. Since PR[#XXX](https://github.com/gridap/GridapDistributed.jl/pull/XXX). + ## [0.4.10] - 2025-09-29 ### Added diff --git a/src/FESpaces.jl b/src/FESpaces.jl index 0a471b69..5871692b 100644 --- a/src/FESpaces.jl +++ b/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) diff --git a/test/FESpacesTests.jl b/test/FESpacesTests.jl index bf1e78e6..76cf57cb 100644 --- a/test/FESpacesTests.jl +++ b/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) From 9121c89cbeb721ca2a7500cc17cc7abc830a6e43 Mon Sep 17 00:00:00 2001 From: "Alberto F. Martin" <38347633+amartinhuertas@users.noreply.github.com> Date: Sun, 2 Nov 2025 17:23:44 +1100 Subject: [PATCH 2/2] Update NEWS.md Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index b1eb72e9..86a11fff 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,7 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -- New overloads for the `TrialFESpace` constructor where the data to be imposed is passed as a `DistributedCellField`. Since PR[#XXX](https://github.com/gridap/GridapDistributed.jl/pull/XXX). +- 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