diff --git a/NEWS.md b/NEWS.md index c38d12a72..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 @@ -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 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" 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 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..99396130c 100644 --- a/src/FESpaces/ZeroMeanFESpaces.jl +++ b/src/FESpaces/ZeroMeanFESpaces.jl @@ -26,16 +26,29 @@ 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, 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,7 +70,16 @@ 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 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 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/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..46b298c0e 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 @@ -31,15 +32,79 @@ 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) +# 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 -mean1 = sum(∫(uh)*dΩ) +# 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 -tol = 1.0e-10 -@test abs(mean1) < tol +# Check Stokes/Navier-Stokes problem properties -V = FESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:L2,constraint=:zeromean) -@test isa(V,ZeroMeanFESpace) +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(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