diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 25954a8d..a14c19f2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -25,7 +25,7 @@ jobs: arch: ${{ matrix.arch }} - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl poisson.jl validation.jl elasticity.jl + - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl poisson.jl validation.jl elasticity.jl lagrange_multipliers.jl tutorials_set_2: name: Tutorials2 ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} runs-on: ${{ matrix.os }} @@ -46,7 +46,7 @@ jobs: arch: ${{ matrix.arch }} - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl p_laplacian.jl hyperelasticity.jl dg_discretization.jl fsi_tutorial.jl + - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl p_laplacian.jl hyperelasticity.jl dg_discretization.jl fsi_tutorial.jl poisson_amr.jl tutorials_set_3: name: Tutorials3 ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} runs-on: ${{ matrix.os }} diff --git a/Project.toml b/Project.toml index 374b111b..506480e2 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.17.0" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e" @@ -33,6 +34,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] +DataStructures = "0.18.22" Gridap = "0.18" GridapDistributed = "0.4" GridapGmsh = "0.7" diff --git a/deps/build.jl b/deps/build.jl index 1470ef80..777cbf82 100644 --- a/deps/build.jl +++ b/deps/build.jl @@ -26,6 +26,7 @@ files = [ "Topology optimization"=>"TopOptEMFocus.jl", "Poisson on unfitted meshes"=>"poisson_unfitted.jl", "Poisson with AMR"=>"poisson_amr.jl", + "Lagrange multipliers" => "lagrange_multipliers.jl", "Low-level API - Poisson equation"=>"poisson_dev_fe.jl", "Low-level API - Geometry" => "geometry_dev.jl", ] diff --git a/src/lagrange_multipliers.jl b/src/lagrange_multipliers.jl new file mode 100644 index 00000000..a9af1031 --- /dev/null +++ b/src/lagrange_multipliers.jl @@ -0,0 +1,127 @@ +# In this tutorial, we will learn +# +# - How to enforce constraints using Lagrange multipliers +# - How to work with `ConstantFESpace` +# +# ## Problem statement +# +# In this tutorial, we solve the Poisson equation with pure Neumann boundary conditions. +# This problem is well-known to be singular since the solution is defined up to a constant, +# which we have to fix to obtain a unique solution. +# Here, we will use a Lagrange multiplier to enforce that the mean value of the solution +# equals a given constant. +# +# The problem reads: find $u$ and $λ$ such that +# +# ```math +# \left\lbrace +# \begin{aligned} +# -\Delta u = f \ &\text{in} \ \Omega,\\ +# \nabla u\cdot n = g \ &\text{on}\ \Gamma,\\ +# \int_{\Omega} u \ {\rm d}\Omega = \bar{u},\\ +# \end{aligned} +# \right. +# ``` +# +# where $\Omega$ is our domain, $\Gamma$ is its boundary, $n$ is the outward unit normal vector, +# and $\bar{u}$ is a given constant that fixes the mean value of the solution. +# +# ## Numerical scheme +# +# The weak form of this problem using Lagrange multipliers reads: +# find $(u,λ) \in V \times \Lambda$ such that +# +# ```math +# \begin{aligned} +# \int_{\Omega} \nabla u \cdot \nabla v \ {\rm d}\Omega + +# \int_{\Omega} λv \ {\rm d}\Omega + +# \int_{\Omega} uμ \ {\rm d}\Omega = +# \int_{\Omega} fv \ {\rm d}\Omega + +# \int_{\Gamma} v(g\cdot n) \ {\rm d}\Gamma + +# \int_{\Omega} μ\bar{u} \ {\rm d}\Omega +# \end{aligned} +# ``` +# +# for all $(v,μ) \in V \times \Lambda$, where $V = H^1(\Omega)$ and $\Lambda = \mathbb{R}$. +# +# ## Implementation +# +# First, we load the Gridap package and define the exact solution that we will use to +# manufacture the source term and boundary condition: + +using Gridap + +u_exact(x) = sin(x[1]) * cos(x[2]) + +# Now we can create a simple Cartesian mesh of the unit square: + +model = CartesianDiscreteModel((0,1,0,1),(8,8)) + +# We will use first order Lagrangian finite elements for the primal variable u. + +order = 1 +reffe = ReferenceFE(lagrangian, Float64, order) +V = FESpace(model, reffe) + +# For the Lagrange multiplier λ, we need a space of constant functions, since λ ∈ ℝ. +# In Gridap, we can create such a space using `ConstantFESpace`: + +Λ = ConstantFESpace(model) + +# Conceptually, a `ConstantFESpace` is a space defined on the whole domain with a +# single degree of freedom, which is what we need for the Lagrange multiplier λ. +# We finally bundle both spaces into a multi-field space: + +X = MultiFieldFESpace([V, Λ]) + +# ## Integration +# +# We need to create the triangulation and measures for both domain and boundary +# integration: + +Ω = Triangulation(model) +Γ = BoundaryTriangulation(model) +dΩ = Measure(Ω, 2*order) +dΓ = Measure(Γ, 2*order) + +# Next, we manufacture the source term f and Neumann boundary condition g +# from the exact solution. We also compute the mean value ū that we want +# to enforce: + +f(x) = -Δ(u_exact)(x) +g(x) = ∇(u_exact)(x) +ū = sum(∫(u_exact)dΩ) +nΓ = get_normal_vector(Γ) + +# ## Weak Form +# +# We can now define the bilinear and linear forms of our problem. +# Note how the forms take tuples as arguments, representing the +# multi-field nature of our solution: + +a((u,λ),(v,μ)) = ∫(∇(u)⋅∇(v) + λ*v + u*μ)dΩ +l((v,μ)) = ∫(f*v + μ*ū)dΩ + ∫(v*(g⋅nΓ))*dΓ + +# ## Solution +# +# We can now create the FE operator and solve the system: + +op = AffineFEOperator(a, l, X, X) +uh, λh = solve(op) + +# Note how we get two values from solve: the primal solution uh and +# the Lagrange multiplier λh. Finally, we compute the L2 error and +# verify that the mean value constraint is satisfied: + +eh = uh - u_exact +l2_error = sqrt(sum(∫(eh⋅eh)*dΩ)) +ūh = sum(∫(uh)*dΩ) + +# The L2 error should be small (of order h²) and ūh should be very close to ū, +# showing that both the equation and the constraint are well satisfied. + +# ## Visualization +# +# We can visualize the solution and error by writing them to a VTK file: + +writevtk(Ω, "results", cellfields=["uh"=>uh, "error"=>eh]) diff --git a/src/poisson_amr.jl b/src/poisson_amr.jl index 6d0deb16..564aab24 100644 --- a/src/poisson_amr.jl +++ b/src/poisson_amr.jl @@ -79,7 +79,11 @@ function LShapedModel(n) cell_coords = map(mean,get_cell_coordinates(model)) l_shape_filter(x) = (x[1] < 0.5) || (x[2] < 0.5) mask = map(l_shape_filter,cell_coords) - return simplexify(DiscreteModelPortion(model,mask)) + model = simplexify(DiscreteModelPortion(model,mask)) + + grid = get_grid(model) + topo = get_grid_topology(model) + return UnstructuredDiscreteModel(grid, topo, FaceLabeling(topo)) end # Define the L2 norm for error estimation. @@ -179,8 +183,8 @@ for i in 1:nsteps ) println("Error: $error, Error η: $(sum(η))") - last_error = error - model = fmodel + global last_error = error + global model = fmodel end # The final mesh gives the following result: diff --git a/test/runtests.jl b/test/runtests.jl index 2d9b398a..7456f936 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,23 +17,17 @@ else end for (title,filename) in files - # Create temporal modules to isolate and protect test scopes - tmpdir=mktempdir(;cleanup=true) - filename_wo_extension=split(filename,".")[1] - tmpmod = filename_wo_extension - tmpfile = joinpath(tmpdir,tmpmod) - isfile(tmpfile) && error("File $tmpfile already exists!") - testpath = joinpath(@__DIR__,"../src", filename) - open(tmpfile,"w") do f - println(f, "# This file is automatically generated") - println(f, "# Do not edit") - println(f) - println(f, "module $tmpmod include(\"$testpath\") end") - end - @time @testset "$title" begin include(tmpfile) end + # Create temporal modules to isolate and protect test scopes + tmpdir = mktempdir(;cleanup=true) + tmpmod = split(filename,".")[1] + tmpfile = joinpath(tmpdir,tmpmod) + isfile(tmpfile) && error("File $tmpfile already exists!") + testpath = joinpath(@__DIR__,"../src", filename) + open(tmpfile,"w") do f + println(f, "# This file is automatically generated") + println(f, "# Do not edit") + println(f) + println(f, "module $tmpmod include(\"$testpath\") end") + end + @time @testset "$title" begin include(tmpfile) end end - -# module fsi_tutorial -# using Test -# @time @testset "fsi_tutorial" begin include("../src/fsi_tutorial.jl") end -# end # module