{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Tutorial 21: Poisson with AMR"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"In this tutorial, we will learn:\n",
"\n",
" - How to use adaptive mesh refinement (AMR) with Gridap\n",
" - How to set up a Poisson problem on an L-shaped domain\n",
" - How to implement error estimation and marking strategies\n",
" - How to visualize the AMR process and results\n",
"\n",
"## Problem Overview\n",
"\n",
"We will solve the Poisson equation on an L-shaped domain using adaptive mesh refinement.\n",
"The L-shaped domain is known to introduce a singularity at its reentrant corner,\n",
"making it an excellent test case for AMR. The problem is:\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"-\\Delta u &= f &&\\text{in } \\Omega \\\\\n",
"u &= g &&\\text{on } \\partial\\Omega\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"where Ω is the L-shaped domain, and we choose an exact solution with a singularity\n",
"to demonstrate the effectiveness of AMR.\n",
"\n",
"## Error Estimation and AMR Process\n",
"\n",
"We use a residual-based a posteriori error estimator η. For each element K in the mesh,\n",
"the local error indicator ηK is computed as:\n",
"\n",
"$$\n",
"\\eta_K^2 = h_K^2\\|f + \\Delta u_h\\|_{L^2(K)}^2 + h_K\\|\\lbrack\\!\\lbrack \\nabla u_h \\cdot n \\rbrack\\!\\rbrack \\|_{L^2(\\partial K)}^2\n",
"$$\n",
"\n",
"where:\n",
"- `h_K` is the diameter of element K\n",
"- `u_h` is the computed finite element solution\n",
"- The first term measures the element residual\n",
"- The second term measures the jump in the normal derivative across element boundaries\n",
"\n",
"The AMR process follows these steps in each iteration:\n",
"\n",
"1. **Solve**: Compute the finite element solution u_h on the current mesh\n",
"2. **Estimate**: Calculate error indicators ηK for each element\n",
"3. **Mark**: Use Dörfler marking to select elements for refinement\n",
" - Sort elements by error indicator\n",
" - Mark elements containing a fixed fraction (here 80%) of total error\n",
"4. **Refine**: Refine selected elements to obtain a new mesh. In this example,\n",
" we will be using the newest vertex bisection (NVB) method to keep the mesh\n",
" conforming (without any hanging nodes).\n",
"\n",
"This adaptive loop continues until either:\n",
"- A maximum number of iterations is reached\n",
"- The estimated error falls below a threshold\n",
"- The solution achieves desired accuracy\n",
"\n",
"The process automatically concentrates mesh refinement in regions of high error,\n",
"particularly around the reentrant corner where the solution has reduced regularity.\n",
"This results in better accuracy per degree of freedom compared to uniform refinement.\n",
"\n",
"## Required Packages"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using Gridap, Gridap.Geometry, Gridap.Adaptivity\n",
"using DataStructures"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"## Problem Setup\n",
"\n",
"We define an exact solution that contains a singularity at the corner (0.5, 0.5)\n",
"of the L-shaped domain. This singularity will demonstrate how AMR automatically\n",
"refines the mesh in regions of high error."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"ϵ = 1e-2\n",
"r(x) = ((x[1]-0.5)^2 + (x[2]-0.5)^2)^(1/2)\n",
"u_exact(x) = 1.0 / (ϵ + r(x))"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"Create an L-shaped domain by removing a quadrant from a unit square.\n",
"The domain is [0,1]² \\ [0.5,1]×[0.5,1]"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function LShapedModel(n)\n",
" model = CartesianDiscreteModel((0,1,0,1),(n,n))\n",
" cell_coords = map(mean,get_cell_coordinates(model))\n",
" l_shape_filter(x) = (x[1] < 0.5) || (x[2] < 0.5)\n",
" mask = map(l_shape_filter,cell_coords)\n",
" return simplexify(DiscreteModelPortion(model,mask))\n",
"end"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"Define the L2 norm for error estimation.\n",
"These will be used to compute both local and global error measures."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"l2_norm(he,xh,dΩ) = ∫(he*(xh*xh))*dΩ\n",
"l2_norm(xh,dΩ) = ∫(xh*xh)*dΩ"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"## AMR Step Function\n",
"\n",
"The `amr_step` function performs a single step of the adaptive mesh refinement process:\n",
"1. Solves the Poisson problem on the current mesh\n",
"2. Estimates the error using residual-based error indicators\n",
"3. Marks cells for refinement using Dörfler marking\n",
"4. Refines the mesh using newest vertex bisection (NVB)"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function amr_step(model,u_exact;order=1)\n",
" \"Create FE spaces with Dirichlet boundary conditions on all boundaries\"\n",
" reffe = ReferenceFE(lagrangian,Float64,order)\n",
" V = TestFESpace(model,reffe;dirichlet_tags=[\"boundary\"])\n",
" U = TrialFESpace(V,u_exact)\n",
"\n",
" \"Setup integration measures\"\n",
" Ω = Triangulation(model)\n",
" Γ = Boundary(model)\n",
" Λ = Skeleton(model)\n",
"\n",
" dΩ = Measure(Ω,4*order)\n",
" dΓ = Measure(Γ,2*order)\n",
" dΛ = Measure(Λ,2*order)\n",
"\n",
" \"Compute cell sizes for error estimation\"\n",
" hK = CellField(sqrt.(collect(get_array(∫(1)dΩ))),Ω)\n",
"\n",
" \"Get normal vectors for boundary and interface terms\"\n",
" nΓ = get_normal_vector(Γ)\n",
" nΛ = get_normal_vector(Λ)\n",
"\n",
" \"Define the weak form\"\n",
" ∇u(x) = ∇(u_exact)(x)\n",
" f(x) = -Δ(u_exact)(x)\n",
" a(u,v) = ∫(∇(u)⋅∇(v))dΩ\n",
" l(v) = ∫(f*v)dΩ\n",
"\n",
" \"Define the residual error estimator\n",
" It includes volume residual, boundary jump, and interface jump terms\"\n",
" ηh(u) = l2_norm(hK*(f + Δ(u)),dΩ) + # Volume residual\n",
" l2_norm(hK*(∇(u) - ∇u)⋅nΓ,dΓ) + # Boundary jump\n",
" l2_norm(jump(hK*∇(u)⋅nΛ),dΛ) # Interface jump\n",
"\n",
" \"Solve the FE problem\"\n",
" op = AffineFEOperator(a,l,U,V)\n",
" uh = solve(op)\n",
"\n",
" \"Compute error indicators\"\n",
" η = estimate(ηh,uh)\n",
"\n",
" \"Mark cells for refinement using Dörfler marking\n",
" This strategy marks cells containing a fixed fraction (0.9) of the total error\"\n",
" m = DorflerMarking(0.9)\n",
" I = Adaptivity.mark(m,η)\n",
"\n",
" \"Refine the mesh using newest vertex bisection\"\n",
" method = Adaptivity.NVBRefinement(model)\n",
" amodel = refine(method,model;cells_to_refine=I)\n",
" fmodel = Adaptivity.get_model(amodel)\n",
"\n",
" \"Compute the global error for convergence testing\"\n",
" error = sum(l2_norm(uh - u_exact,dΩ))\n",
" return fmodel, uh, η, I, error\n",
"end"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"## Main AMR Loop\n",
"\n",
"We perform multiple AMR steps, refining the mesh iteratively and solving\n",
"the problem on each refined mesh. This demonstrates how the error decreases\n",
"as the mesh is adaptively refined in regions of high error."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"nsteps = 5\n",
"order = 1\n",
"model = LShapedModel(10)\n",
"\n",
"last_error = Inf\n",
"for i in 1:nsteps\n",
" fmodel, uh, η, I, error = amr_step(model,u_exact;order)\n",
"\n",
" is_refined = map(i -> ifelse(i ∈ I, 1, -1), 1:num_cells(model))\n",
"\n",
" Ω = Triangulation(model)\n",
" writevtk(\n",
" Ω,\"model_$(i-1)\",append=false,\n",
" cellfields = [\n",
" \"uh\" => uh, # Computed solution\n",
" \"η\" => CellField(η,Ω), # Error indicators\n",
" \"is_refined\" => CellField(is_refined,Ω), # Refinement markers\n",
" \"u_exact\" => CellField(u_exact,Ω), # Exact solution\n",
" ],\n",
" )\n",
"\n",
" println(\"Error: $error, Error η: $(sum(η))\")\n",
" last_error = error\n",
" model = fmodel\n",
"end"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"The final mesh gives the following result:\n",
"\n",
"\n",
"## Conclusion\n",
"\n",
"In this tutorial, we have demonstrated how to:\n",
"1. Implement adaptive mesh refinement for the Poisson equation\n",
"2. Use residual-based error estimation to identify regions needing refinement\n",
"3. Apply Dörfler marking to select cells for refinement\n",
"4. Visualize the AMR process and solution convergence\n",
"\n",
"The results show how AMR automatically refines the mesh near the singularity,\n",
"leading to more efficient and accurate solutions compared to uniform refinement."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"---\n",
"\n",
"*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
],
"metadata": {}
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.10.8"
},
"kernelspec": {
"name": "julia-1.10",
"display_name": "Julia 1.10.8",
"language": "julia"
}
},
"nbformat": 4
}