Merge pull request #1020 from gridap/distributed-zero-mean

ZeroMeanFESpace bugfixes

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
 

Project.toml

@@ -1,7 +1,7 @@
 name = "Gridap"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 authors = ["Santiago Badia <[email protected]>", "Francesc Verdugo <[email protected]>", "Alberto F. Martin <[email protected]>"]
-version = "0.18.3"
+version = "0.18.4"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

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

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)

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

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)

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

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