Add MeanValueRestriction

Author: pjaap

examples/Example252_NSEPlanarLatticeFlow.jl

@@ -83,7 +83,14 @@ function exact_error!(result, u, qpinfo)
     return result .= result .^ 2
 end
 
-function main(; μ = 0.001, nrefs = 5, reconstruct = true, Plotter = nothing, kwargs...)
+function main(;
+        μ = 0.001,
+        nrefs = 5,
+        reconstruct = true,
+        Plotter = nothing,
+        use_LM_restrictions = true,
+        kwargs...
+    )
 
     ## problem description
     PD = ProblemDescription()
@@ -96,6 +103,11 @@ function main(; μ = 0.001, nrefs = 5, reconstruct = true, Plotter = nothing, kw
     assign_operator!(PD, NonlinearOperator(kernel_nonlinear!, [id_u, grad(u), id(p)]; params = [μ], kwargs...))
     assign_operator!(PD, LinearOperator(f!(μ), [id_u]; kwargs...))
     assign_operator!(PD, InterpolateBoundaryData(u, u!; regions = 1:4))
+    if use_LM_restrictions
+        assign_restriction!(PD, MeanValueRestriction(p; value = 0))
+    else
+        assign_operator!(PD, FixDofs(p; dofs = [1]))
+    end
 
     ## grid
     xgrid = uniform_refine(grid_unitsquare(Triangle2D), nrefs)
@@ -106,10 +118,12 @@ function main(; μ = 0.001, nrefs = 5, reconstruct = true, Plotter = nothing, kw
     ## solve
     sol = solve(PD, FES; kwargs...)
 
-    ## move integral mean of pressure
-    pintegrate = ItemIntegrator([id(p)])
-    pmean = sum(evaluate(pintegrate, sol)) / sum(xgrid[CellVolumes])
-    view(sol[p]) .-= pmean
+    if !use_LM_restrictions
+        ## move integral mean of pressure
+        pintegrate = ItemIntegrator([id(p)])
+        pmean = sum(evaluate(pintegrate, sol)) / sum(xgrid[CellVolumes])
+        view(sol[p]) .-= pmean
+    end
 
     ## error calculation
     ErrorIntegratorExact = ItemIntegrator(exact_error!, [id(u), id(p)]; quadorder = 4, params = [μ], kwargs...)
@@ -127,7 +141,9 @@ end
 
 generateplots = ExtendableFEM.default_generateplots(Example252_NSEPlanarLatticeFlow, "example252.png") #hide
 function runtests() #hide
-    L2errorU, plt = main(; nrefs = 4) #hide
+    L2errorU, _ = main(; use_LM_restrictions = false, nrefs = 4) #hide
+    @test L2errorU ≈ 0.010616923333947861 #hide
+    L2errorU, _ = main(; use_LM_restrictions = true, nrefs = 4) #hide
     @test L2errorU ≈ 0.010616923333947861 #hide
     return nothing #hide
 end #hide

src/ExtendableFEM.jl

@@ -128,6 +128,8 @@ include("common_operators/reduction_operator.jl")
 include("restrictions.jl")
 include("common_restrictions/coupled_dofs_restriction.jl")
 export CoupledDofsRestriction
+include("common_restrictions/mean_value_restriction.jl")
+export MeanValueRestriction
 
 include("problemdescription.jl")
 export ProblemDescription

src/common_restrictions/coupled_dofs_restriction.jl

@@ -20,7 +20,7 @@ function CoupledDofsRestriction(matrix::AbstractMatrix)
 end
 
 
-function assemble!(R::CoupledDofsRestriction, SC; kwargs...)
+function assemble!(R::CoupledDofsRestriction, sol, SC; kwargs...)
 
     # extract all col indices
     _, J, _ = findnz(R.coupling_matrix)

src/common_restrictions/mean_value_restriction.jl

@@ -0,0 +1,36 @@
+struct MeanValueRestriction{T, Tv} <: AbstractRestriction
+    unknown::Unknown
+    value::T
+    linear_operator
+    parameters::Dict{Symbol, Any}
+end
+
+
+"""
+$(TYPEDSIGNATURES)
+
+Restrict the mean value of a scalar unknown to a certain value.
+
+An additional operator can be applied to the unknown before the mean value is restricted.
+"""
+function MeanValueRestriction(u::Unknown; value::T = 0, operator = Identity, Tv = Float64) where {T}
+    linear_operator = LinearOperator([(u, operator)]; store = true, Tv)
+    return MeanValueRestriction{T, Tv}(u, value, linear_operator, Dict{Symbol, Any}(:name => "MeanValueRestriction"))
+end
+
+
+function assemble!(R::MeanValueRestriction{T, Tv}, sol, SC; kwargs...) where {T, Tv}
+
+    n = length(SC.b.entries)
+
+    # type and shape of the rhs
+    b = deepcopy(SC.b)
+    fill!(b.entries, 0)
+
+    assemble!(nothing, b, sol, R.linear_operator, SC; assemble_rhs = true, kwargs...)
+
+    R.parameters[:matrix] = reshape(b.entries, n, 1)
+    R.parameters[:rhs] = Tv[R.value]
+
+    return nothing
+end

src/solvers.jl

@@ -227,7 +227,7 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
         ## assemble restrctions
         if !SC.parameters[:initialized]
             for restriction in PD.restrictions
-                @timeit timer "$(restriction.parameters[:name])" assemble!(restriction, SC; kwargs...)
+                @timeit timer "$(restriction.parameters[:name])" assemble!(restriction, sol, SC; kwargs...)
             end
         end
     end
@@ -287,9 +287,9 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
                 end
 
                 b_block = BlockVector(zeros(Tv, total_size), block_sizes)
-                b_block[Block(1)] = b_unrestricted
 
-                u_unrestricted = linsolve.u
+                b_block[Block(1)] = b_unrestricted
+                u_unrestricted = @views linsolve.u[1:block_sizes[1]]
                 u_block = BlockVector(zeros(Tv, total_size), block_sizes)
                 u_block[Block(1)] = u_unrestricted
 
@@ -333,7 +333,8 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
         # Solve linear system
         push!(stats[:matrix_nnz], nnz(linsolve.A))
         @timeit timer "solve! call" begin
-            blocked_Δx = LinearSolve.solve!(linsolve)
+            blocked_result = LinearSolve.solve!(linsolve)
+            blocked_Δx = blocked_result.u
         end
 
         # extract the solution / dismiss the lagrange multipliers