From 904c52a74a31cc6dc7e4776b5a76864074931b12 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 7 Aug 2024 17:40:55 +1000 Subject: [PATCH 1/8] Bugfix: FESpaceWithConstantFixed was not modifying free values --- src/FESpaces/FESpacesWithConstantFixed.jl | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/FESpaces/FESpacesWithConstantFixed.jl b/src/FESpaces/FESpacesWithConstantFixed.jl index a76f54260..95156ee8f 100644 --- a/src/FESpaces/FESpacesWithConstantFixed.jl +++ b/src/FESpaces/FESpacesWithConstantFixed.jl @@ -45,15 +45,12 @@ function get_cell_dof_ids(f::FESpaceWithConstantFixed{DoNotFixConstant}) end get_dirichlet_dof_ids(f::FESpaceWithConstantFixed{FixConstant}) = Base.OneTo(1) - get_dirichlet_dof_ids(f::FESpaceWithConstantFixed{DoNotFixConstant}) = Base.OneTo(0) num_dirichlet_tags(f::FESpaceWithConstantFixed{FixConstant}) = 1 - num_dirichlet_tags(f::FESpaceWithConstantFixed{DoNotFixConstant}) = 0 get_dirichlet_dof_tag(f::FESpaceWithConstantFixed{FixConstant}) = Int8[1,] - get_dirichlet_dof_tag(f::FESpaceWithConstantFixed{DoNotFixConstant}) = Int8[] function scatter_free_and_dirichlet_values(f::FESpaceWithConstantFixed{FixConstant},fv,dv) @@ -81,15 +78,16 @@ function gather_free_and_dirichlet_values(f::FESpaceWithConstantFixed{DoNotFixCo end function gather_free_and_dirichlet_values!(fv,dv,f::FESpaceWithConstantFixed{FixConstant},cv) - _fv, _dv = gather_free_and_dirichlet_values(f.space,cv) - @assert length(_dv) == 0 - fv .= VectorWithEntryRemoved(_fv,f.dof_to_fix) - dv[1] = _fv[f.dof_to_fix] + @assert length(dv) == 1 + _dv = similar(dv,eltype(dv),0) + _fv = VectorWithEntryInserted(fv,f.dof_to_fix,zero(eltype(fv))) + gather_free_and_dirichlet_values!(_fv,_dv,f.space,cv) + dv[1] = _fv.value (fv, dv) end function gather_free_and_dirichlet_values!(fv,dv,f::FESpaceWithConstantFixed{DoNotFixConstant},cv) - gather_free_and_dirichlet_values(f.space,cv) + gather_free_and_dirichlet_values!(fv,dv,f.space,cv) end function TrialFESpace(f::FESpaceWithConstantFixed{CA}) where CA From 8c292162ed76e2fdb27b02aa4a432f16aac00f81 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 7 Aug 2024 18:19:29 +1000 Subject: [PATCH 2/8] Bugfix: ZeroMeanFESpace was not interpolating correctly --- src/FESpaces/SingleFieldFESpaces.jl | 6 +++--- src/FESpaces/ZeroMeanFESpaces.jl | 9 +++++++++ 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/src/FESpaces/SingleFieldFESpaces.jl b/src/FESpaces/SingleFieldFESpaces.jl index 674c9e209..2e3538608 100644 --- a/src/FESpaces/SingleFieldFESpaces.jl +++ b/src/FESpaces/SingleFieldFESpaces.jl @@ -184,9 +184,9 @@ end """ """ function interpolate!(object, free_values,fs::SingleFieldFESpace) - cell_vals = _cell_vals(fs,object) - gather_free_values!(free_values,fs,cell_vals) - FEFunction(fs,free_values) + 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) diff --git a/src/FESpaces/ZeroMeanFESpaces.jl b/src/FESpaces/ZeroMeanFESpaces.jl index f1bc563b9..2549adead 100644 --- a/src/FESpaces/ZeroMeanFESpaces.jl +++ b/src/FESpaces/ZeroMeanFESpaces.jl @@ -60,6 +60,15 @@ function _compute_new_fixedval(fv,dv,vol_i,vol,fixed_dof) c end +# This is required, otherwise we end up calling `FEFunction` with a fixed value of zero, +# which does not properly interpolate the function provided. +# With this change, we modify are interpolating in the unconstrained space and then +# substracting the mean. +function interpolate!(object,free_values,fs::ZeroMeanFESpace) + dirichlet_values = zero_dirichlet_values(fs) + interpolate_everywhere!(object,free_values,dirichlet_values,fs) +end + # Delegated functions get_triangulation(f::ZeroMeanFESpace) = get_triangulation(f.space) From bd8c530a780228035bdad3712d1bc40ef8355ba6 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 7 Aug 2024 23:39:11 +1000 Subject: [PATCH 3/8] Added tests --- .../FESpacesWithConstantFixedTests.jl | 10 ++++++++-- test/FESpacesTests/ZeroMeanFESpacesTests.jl | 20 ++++++++++++------- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/test/FESpacesTests/FESpacesWithConstantFixedTests.jl b/test/FESpacesTests/FESpacesWithConstantFixedTests.jl index 68c5cb999..da22ce92b 100644 --- a/test/FESpacesTests/FESpacesWithConstantFixedTests.jl +++ b/test/FESpacesTests/FESpacesWithConstantFixedTests.jl @@ -23,8 +23,14 @@ test_single_field_fe_space(V0) uh0 = interpolate(V0) do x sin(4*pi*(x[1]+x[2]^2)) + 3 end -using Gridap.Visualization -#writevtk(trian,"trian",nsubcells=20,cellfields=["uh0"=>uh0]) +V1 = FESpaceWithConstantFixed(V,false,rand(1:num_free_dofs(V))) +test_single_field_fe_space(V1) + +@test Gridap.FESpaces.ConstantApproach(V1) == Gridap.FESpaces.DoNotFixConstant() + +uh1 = interpolate(V1) do x + sin(4*pi*(x[1]+x[2]^2)) + 3 +end end # module diff --git a/test/FESpacesTests/ZeroMeanFESpacesTests.jl b/test/FESpacesTests/ZeroMeanFESpacesTests.jl index ed78ea91e..4f98ecc5f 100644 --- a/test/FESpacesTests/ZeroMeanFESpacesTests.jl +++ b/test/FESpacesTests/ZeroMeanFESpacesTests.jl @@ -31,13 +31,19 @@ U = TrialFESpace(V) test_single_field_fe_space(U,cellmatvec,cellmat,cellvec,trian) @test isa(U,ZeroMeanFESpace) -fun(x) = sin(4*pi*(x[1]+x[2]^2)) + 3 -uh = interpolate(fun, U) - -mean1 = sum(∫(uh)*dΩ) - -tol = 1.0e-10 -@test abs(mean1) < tol +# Interpolate a function with non-zero mean +f(x) = sin(4*pi*(x[1]+x[2]^2)) + 3 +uh = interpolate(f, U) +@test abs(sum(∫(uh)*dΩ)) < 1.0e-10 + +# Interpolate a function with zero mean +ĝ(x) = x[1]^2 + x[2] +g_mean = sum(∫(ĝ)*dΩ)/sum(∫(1)*dΩ) +g(x) = ĝ(x) - g_mean +vh = interpolate(g, U) +eh = vh - g +@test abs(sum(∫(vh)*dΩ)) < 1.0e-10 +@test sum(∫(eh*eh)*dΩ) < 1.0e-10 V = FESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:L2,constraint=:zeromean) @test isa(V,ZeroMeanFESpace) From 629e9cdab369516517dbcd8406204dd3a5bc8da0 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 8 Aug 2024 07:49:29 +1000 Subject: [PATCH 4/8] Updated NEWS.md --- NEWS.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NEWS.md b/NEWS.md index c38d12a72..962891492 100644 --- a/NEWS.md +++ b/NEWS.md @@ -14,6 +14,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed - Passing `kwargs` from `refine` to `simplexify` functions in Adaptivity. Since PR[#1015](https://github.com/gridap/Gridap.jl/pull/1015). +- Fixed `interpolate` for `ZeroMeanFESpace`. Since PR[#1020](https://github.com/gridap/Gridap.jl/pull/1020). +- Fixed `gather_free_and_dirichlet_values!` for `FESpaceWithConstantFixed`. Since PR[#1020](https://github.com/gridap/Gridap.jl/pull/1020). ## [0.18.3] - 2024-07-11 From 86a02fc9da9bc43baba8553a71b505622a828074 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 8 Aug 2024 13:05:39 +1000 Subject: [PATCH 5/8] Added more tests for ZeroMeanFESpaces --- src/FESpaces/ZeroMeanFESpaces.jl | 7 ++-- test/FESpacesTests/ZeroMeanFESpacesTests.jl | 42 ++++++++++++++++++++- 2 files changed, 44 insertions(+), 5 deletions(-) diff --git a/src/FESpaces/ZeroMeanFESpaces.jl b/src/FESpaces/ZeroMeanFESpaces.jl index 2549adead..12e46172b 100644 --- a/src/FESpaces/ZeroMeanFESpaces.jl +++ b/src/FESpaces/ZeroMeanFESpaces.jl @@ -35,7 +35,8 @@ function FEFunction( dirichlet_values, f.vol_i, f.vol, - f.space.dof_to_fix) + f.space.dof_to_fix + ) fv = lazy_map(+,free_values,Fill(c,length(free_values))) dv = dirichlet_values .+ c FEFunction(f.space,fv,dv) @@ -57,12 +58,12 @@ function _compute_new_fixedval(fv,dv,vol_i,vol,fixed_dof) c += fv[i-1]*vol_i[i] end c = -c/vol - c + return c end # This is required, otherwise we end up calling `FEFunction` with a fixed value of zero, # which does not properly interpolate the function provided. -# With this change, we modify are interpolating in the unconstrained space and then +# With this change, we are interpolating in the unconstrained space and then # substracting the mean. function interpolate!(object,free_values,fs::ZeroMeanFESpace) dirichlet_values = zero_dirichlet_values(fs) diff --git a/test/FESpacesTests/ZeroMeanFESpacesTests.jl b/test/FESpacesTests/ZeroMeanFESpacesTests.jl index 4f98ecc5f..e37a3266f 100644 --- a/test/FESpacesTests/ZeroMeanFESpacesTests.jl +++ b/test/FESpacesTests/ZeroMeanFESpacesTests.jl @@ -1,6 +1,7 @@ module ZeroMeanFESpacesTests using Test +using Gridap using Gridap.Arrays using Gridap.Geometry using Gridap.Fields @@ -45,7 +46,44 @@ eh = vh - g @test abs(sum(∫(vh)*dΩ)) < 1.0e-10 @test sum(∫(eh*eh)*dΩ) < 1.0e-10 -V = FESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:L2,constraint=:zeromean) -@test isa(V,ZeroMeanFESpace) +# Check Stokes problem properties + +u_ex(x) = VectorValue(x[2],-x[1]) +p_ex(x) = x[1] + 2*x[2] +p_mean = sum(∫(p_ex)*dΩ) + +order = 2 +reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) + +V = FESpace(model,reffe_u;dirichlet_tags="boundary") +U = TrialFESpace(V,u_ex) +Q = FESpace(model,reffe_p;conformity=:L2) +Q0 = FESpace(model,reffe_p;conformity=:L2,constraint=:zeromean) +@test isa(Q0,ZeroMeanFESpace) + +X = MultiFieldFESpace([U,Q0]) +Y = MultiFieldFESpace([V,Q0]) + +ph_i = interpolate(p_ex, Q0) +@test abs(sum(∫(ph_i)*dΩ)) < 1.0e-10 +@test abs(sum(∫(ph_i - p_ex + p_mean)*dΩ)) < 1.0e-10 + +a((u,p),(v,q)) = ∫(∇(u)⊙∇(v) - (∇⋅v)*p - q*(∇⋅u))dΩ +l((v,q)) = a((u_ex,p_ex),(v,q)) + +op = AffineFEOperator(a,l,X,Y) +uh, ph = solve(op) + +l2_error(u,v) = sqrt(sum(∫((u-v)⋅(u-v))*dΩ)) +@test l2_error(uh,u_ex) < 1.0e-10 +@test abs(sum(∫(ph)*dΩ)) < 1.0e-10 +@test l2_error(ph,ph_i) < 1.0e-10 + +b(u,q) = ∫(q*(∇⋅u))dΩ +B = assemble_vector(q -> b(uh,q),Q) +B0 = assemble_vector(q -> b(u_ex,q),Q0) +@test abs(sum(B)) < 1.0e-10 +@test abs(sum(B0)) < 1.0e-10 end # module From 7fa8398ce5e4be798ec2d0a8d6f08295e393d8de Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 8 Aug 2024 13:06:35 +1000 Subject: [PATCH 6/8] Minor --- test/FESpacesTests/ZeroMeanFESpacesTests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/FESpacesTests/ZeroMeanFESpacesTests.jl b/test/FESpacesTests/ZeroMeanFESpacesTests.jl index e37a3266f..da353c458 100644 --- a/test/FESpacesTests/ZeroMeanFESpacesTests.jl +++ b/test/FESpacesTests/ZeroMeanFESpacesTests.jl @@ -82,7 +82,7 @@ l2_error(u,v) = sqrt(sum(∫((u-v)⋅(u-v))*dΩ)) b(u,q) = ∫(q*(∇⋅u))dΩ B = assemble_vector(q -> b(uh,q),Q) -B0 = assemble_vector(q -> b(u_ex,q),Q0) +B0 = assemble_vector(q -> b(uh,q),Q0) @test abs(sum(B)) < 1.0e-10 @test abs(sum(B0)) < 1.0e-10 From 12832821487284691d609c526178667c485c23fc Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 8 Aug 2024 23:18:58 +1000 Subject: [PATCH 7/8] More tests --- src/FESpaces/ZeroMeanFESpaces.jl | 14 ++++++++++++- test/FESpacesTests/ExtendedFESpacesTests.jl | 2 +- test/FESpacesTests/ZeroMeanFESpacesTests.jl | 23 ++++++++++++++++++++- 3 files changed, 36 insertions(+), 3 deletions(-) diff --git a/src/FESpaces/ZeroMeanFESpaces.jl b/src/FESpaces/ZeroMeanFESpaces.jl index 12e46172b..99396130c 100644 --- a/src/FESpaces/ZeroMeanFESpaces.jl +++ b/src/FESpaces/ZeroMeanFESpaces.jl @@ -26,10 +26,22 @@ function TrialFESpace(f::ZeroMeanFESpace) ZeroMeanFESpace(U,f.vol_i,f.vol) end +# function FEFunction(f::ZeroMeanFESpace,free_values) +# msg = """ +# This function should never be called for ZeroMeanFESpace. Please use +# `FEFunction(f::ZeroMeanFESpace,free_values,dirichlet_values)` instead. +# Reason: +# Without the fixed value (i.e the dirichlet_values), we cannot correctly +# interpolate the free dofs within the space. +# """ +# @unreachable msg +# end +# function FEFunction( f::ZeroMeanFESpace, free_values::AbstractVector, - dirichlet_values::AbstractVector) + dirichlet_values::AbstractVector +) c = _compute_new_fixedval( free_values, dirichlet_values, diff --git a/test/FESpacesTests/ExtendedFESpacesTests.jl b/test/FESpacesTests/ExtendedFESpacesTests.jl index 235fd9d7f..b327656e6 100644 --- a/test/FESpacesTests/ExtendedFESpacesTests.jl +++ b/test/FESpacesTests/ExtendedFESpacesTests.jl @@ -149,7 +149,7 @@ V_in = TestFESpace(Ω_in,ReferenceFE(lagrangian,Float64,2),conformity=:H1) V = TestFESpace(Ω,ReferenceFE(lagrangian,Float64,2),conformity=:H1) vh_in = interpolate(V_in) do x - x[1] + x[1] end vh_in = interpolate(vh_in, V_in) vh = interpolate(vh_in, V) diff --git a/test/FESpacesTests/ZeroMeanFESpacesTests.jl b/test/FESpacesTests/ZeroMeanFESpacesTests.jl index da353c458..46b298c0e 100644 --- a/test/FESpacesTests/ZeroMeanFESpacesTests.jl +++ b/test/FESpacesTests/ZeroMeanFESpacesTests.jl @@ -46,7 +46,7 @@ eh = vh - g @test abs(sum(∫(vh)*dΩ)) < 1.0e-10 @test sum(∫(eh*eh)*dΩ) < 1.0e-10 -# Check Stokes problem properties +# Check Stokes/Navier-Stokes problem properties u_ex(x) = VectorValue(x[2],-x[1]) p_ex(x) = x[1] + 2*x[2] @@ -86,4 +86,25 @@ B0 = assemble_vector(q -> b(uh,q),Q0) @test abs(sum(B)) < 1.0e-10 @test abs(sum(B0)) < 1.0e-10 +uh_ex = interpolate(u_ex, U) +@test l2_error(uh,uh_ex) < 1.0e-10 +conv(u,∇u) = (∇u')⋅u +dconv(du,∇du,u,∇u) = conv(u,∇du)+conv(du,∇u) +c(u,v) = ∫(v⊙(conv∘(u,∇(u))))dΩ +dc(u,du,dv) = ∫(dv⊙(dconv∘(du,∇(du),u,∇(u))))dΩ +jac((u,p),(du,dp),(dv,dq)) = a((du,dp),(dv,dq)) + dc(u,du,dv) +res((u,p),(dv,dq)) = a((u,p),(dv,dq)) + c(u,dv) - a((u_ex,p_ex),(dv,dq)) - c(uh_ex,dv) + +op = FEOperator(res,jac,X,Y) +uh, ph = solve(op) + +@test l2_error(uh,u_ex) < 1.0e-10 +@test abs(sum(∫(ph)*dΩ)) < 1.0e-10 +@test l2_error(ph,ph_i) < 1.0e-10 + +B = assemble_vector(q -> b(uh,q),Q) +B0 = assemble_vector(q -> b(uh,q),Q0) +@test abs(sum(B)) < 1.0e-10 +@test abs(sum(B0)) < 1.0e-10 + end # module From 86c56ded5ae8408d0b368eb9284d72619bae2255 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 9 Aug 2024 10:22:25 +1000 Subject: [PATCH 8/8] Bump version --- NEWS.md | 2 +- Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 962891492..6936b112f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,7 +5,7 @@ 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 +## [0.18.4] - 2024-08-09 ### Changed diff --git a/Project.toml b/Project.toml index 7fd4cb586..69a423369 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Gridap" uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e" authors = ["Santiago Badia ", "Francesc Verdugo ", "Alberto F. Martin "] -version = "0.18.3" +version = "0.18.4" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"