diff --git a/assets/amr/result.png b/assets/amr/result.png
new file mode 100644
index 0000000..15e4b20
Binary files /dev/null and b/assets/amr/result.png differ
diff --git a/assets/geometry/mobius.png b/assets/geometry/mobius.png
new file mode 100644
index 0000000..0fcbcdc
Binary files /dev/null and b/assets/geometry/mobius.png differ
diff --git a/assets/geometry/nodes_fullperiodic.png b/assets/geometry/nodes_fullperiodic.png
new file mode 100644
index 0000000..4922969
Binary files /dev/null and b/assets/geometry/nodes_fullperiodic.png differ
diff --git a/assets/geometry/nodes_halfperiodic.png b/assets/geometry/nodes_halfperiodic.png
new file mode 100644
index 0000000..72410bf
Binary files /dev/null and b/assets/geometry/nodes_halfperiodic.png differ
diff --git a/assets/geometry/nodes_nonperiodic.png b/assets/geometry/nodes_nonperiodic.png
new file mode 100644
index 0000000..2d90167
Binary files /dev/null and b/assets/geometry/nodes_nonperiodic.png differ
diff --git a/deps/build.jl b/deps/build.jl
index 5d1373e..1470ef8 100644
--- a/deps/build.jl
+++ b/deps/build.jl
@@ -18,14 +18,17 @@ files = [
   "Isotropic damage model" => "isotropic_damage.jl",
   "Fluid-Structure Interaction"=>"fsi_tutorial.jl",
   "Electromagnetic scattering in 2D"=>"emscatter.jl",
-  "Low-level API Poisson equation"=>"poisson_dev_fe.jl",
   "On using DrWatson.jl"=>"validation_DrWatson.jl",
   "Interpolation of CellFields"=>"interpolation_fe.jl",
   "Poisson equation on parallel distributed-memory computers"=>"poisson_distributed.jl",
   "Transient Poisson equation"=>"transient_linear.jl",
   "Transient nonlinear equation"=>"transient_nonlinear.jl",
   "Topology optimization"=>"TopOptEMFocus.jl",
-  "Unfitted Poisson"=>"unfitted_poisson.jl"]
+  "Poisson on unfitted meshes"=>"poisson_unfitted.jl",
+  "Poisson with AMR"=>"poisson_amr.jl",
+  "Low-level API - Poisson equation"=>"poisson_dev_fe.jl",
+  "Low-level API - Geometry" => "geometry_dev.jl",
+]
 
 Sys.rm(notebooks_dir;recursive=true,force=true)
 for (i,(title,filename)) in enumerate(files)
diff --git a/src/geometry_dev.jl b/src/geometry_dev.jl
new file mode 100644
index 0000000..100fcf8
--- /dev/null
+++ b/src/geometry_dev.jl
@@ -0,0 +1,391 @@
+#
+# In this tutorial, we will learn
+#
+#    -  How the `DiscreteModel` and its components work.
+#    -  How to extract topological information from a `GridTopology`.
+#    -  How to extract geometrical information from a `Grid`.
+#    -  How periodicity is handled in Gridap, and the difference between nodes and vertices.
+#    -  How to create a periodic model from scratch, use the example of a Mobius strip.
+#
+# ## Required Packages
+
+using Gridap
+using Gridap.Geometry, Gridap.ReferenceFEs, Gridap.Arrays
+using Plots
+
+#
+# ## Table of Contents
+# 1. Utility Functions
+# 2. The DiscreteModel Structure
+# 3. Working with Topology
+# 4. Geometric Mappings
+# 5. High-order Grids
+# 6. Periodicity in Gridap
+#
+# ## 1. Utility Functions
+# We begin by defining helper functions that will be essential throughout this tutorial.
+# These functions help us visualize and work with our mesh structures.
+
+# Convert a CartesianDiscreteModel to an UnstructuredDiscreteModel for more generic handling
+function cartesian_model(args...; kwargs...)
+  UnstructuredDiscreteModel(CartesianDiscreteModel(args...; kwargs...))
+end
+
+# Visualization function to plot nodes with their IDs
+# Input: node_coords - Array of node coordinates
+#        node_ids - Array of corresponding node IDs
+function plot_node_numbering(node_coords, node_ids)
+  x = map(c -> c[1], node_coords)
+  y = map(c -> c[2], node_coords)
+  a = text.(node_ids, halign=:left, valign=:bottom)
+  scatter(x, y, series_annotations = a, legend=false)
+  hline!(unique(x), linestyle=:dash, color=:grey)
+  vline!(unique(y), linestyle=:dash, color=:grey)
+end
+
+# Overloaded method to plot node numbering directly from a model
+# This function extracts the necessary information from the model and calls the base plotting function
+function plot_node_numbering(model)
+  D = num_cell_dims(model)
+  topo = get_grid_topology(model)
+  node_coords = Geometry.get_node_coordinates(model)
+  cell_node_ids = get_cell_node_ids(model)
+  cell_vertex_ids = Geometry.get_faces(topo, D, 0)
+
+  node_to_vertex = zeros(Int, length(node_coords))
+  for (nodes,vertices) in zip(cell_node_ids, cell_vertex_ids)
+    node_to_vertex[nodes] .= vertices
+  end
+  
+  plot_node_numbering(node_coords, node_to_vertex)
+end
+
+#
+# ## 2. The DiscreteModel Structure
+#
+# The `DiscreteModel` in Gridap is a fundamental structure that represents a discretized 
+# computational domain. It consists of three main components, each serving a specific purpose:
+#
+#   - The `GridTopology`: Defines the connectivity of the mesh elements
+#     * Stores how vertices, edges, faces, and cells are connected
+#     * Enables neighbor queries and traversal of the mesh
+#     * Pure topological information, no geometric data
+#
+#   - The `Grid`: Contains the geometric information of the mesh
+#     * Stores coordinates of mesh nodes
+#     * Provides mappings between reference and physical elements
+#     * Handles curved elements and high-order geometries
+#
+#   - The `FaceLabeling`: Manages mesh labels and markers
+#     * Identifies boundary regions
+#     * Tags different material regions
+#     * Essential for applying boundary conditions
+#
+# ### Key Concept: Nodes vs. Vertices
+#
+# One of the most important distinctions in Gridap is between nodes and vertices:
+#
+#   - **Vertices** (Topological entities):
+#     * 0-dimensional entities in the `GridTopology`
+#     * Define the connectivity of the mesh
+#     * Used for neighbor queries and mesh traversal
+#     * Number of vertices depends only on topology
+#
+#   - **Nodes** (Geometric entities):
+#     * Control points stored in the `Grid`
+#     * Define the geometry of elements
+#     * Used for interpolation and mapping
+#     * Number of nodes depends on the geometric order
+#
+# While nodes and vertices often coincide in simple meshes, they differ in two important cases:
+# 1. Periodic meshes: Where multiple nodes may correspond to the same vertex
+# 2. High-order meshes: Where additional nodes are needed to represent curved geometries
+#
+# ## 3. Working with Topology
+#
+# Let's explore how to extract and work with topological information. We'll start by creating
+# a simple 3x3 cartesian mesh:
+
+model = cartesian_model((0,1,0,1),(3,3))
+
+# First, let's get the topology component and count the mesh entities:
+
+topo = get_grid_topology(model)
+
+n_vertices = num_faces(topo,0)  # Number of vertices (0-dimensional entities)
+n_edges = num_faces(topo,1)     # Number of edges (1-dimensional entities)
+n_cells = num_faces(topo,2)     # Number of cells (2-dimensional entities)
+
+# ### Connectivity Queries
+# Gridap provides various methods to query the connectivity between different mesh entities.
+# Here are some common queries:
+
+# Get vertices of each cell (cell → vertex connectivity)
+cell_to_vertices = get_faces(topo,2,0)
+
+# Get edges of each cell (cell → edge connectivity)
+cell_to_edges = get_faces(topo,2,1)
+
+# Get cells adjacent to each edge (edge → cell connectivity)
+edge_to_cells = get_faces(topo,1,2)
+
+# Get vertices of each edge (edge → vertex connectivity)
+edge_to_vertices = get_faces(topo,1,0)
+
+# ### Advanced Connectivity: Finding Cell Neighbors
+#
+# Finding cells that share entities (like vertices or edges) requires more work.
+# A direct query for cell-to-cell connectivity returns an identity map:
+
+cell_to_cells = get_faces(topo,2,2)  # Returns identity table
+
+# To find actual cell neighbors, we need to traverse through lower-dimensional entities.
+# Here's a utility function that builds a face-to-face connectivity graph:
+
+function get_face_to_face_graph(topo,Df)
+  n_faces = num_faces(topo,Df)
+  face_to_vertices = get_faces(topo,Df,0)  # Get vertices of each face
+  vertex_to_faces = get_faces(topo,0,Df)   # Get faces incident to each vertex
+
+  face_to_face = Vector{Vector{Int}}(undef,n_faces)
+  for face in 1:n_faces
+    nbors = Int[]
+    for vertex in face_to_vertices[face]
+      append!(nbors,vertex_to_faces[vertex]) # Add incident faces
+    end
+    face_to_face[face] = filter(!isequal(face),unique(nbors)) # Remove self-reference and duplicates
+  end
+
+  return face_to_face
+end
+
+# Now we can find neighboring cells and edges:
+cell_to_cells = get_face_to_face_graph(topo,2)  # Cells sharing vertices
+edge_to_edges = get_face_to_face_graph(topo,1)  # Edges sharing vertices
+
+#
+# ## 4. Geometric Mappings
+#
+# The geometry of our mesh is defined by mappings from reference elements to physical space.
+# Let's explore how these mappings work in Gridap:
+
+grid = get_grid(model)
+
+# First, we extract the basic geometric information:
+cell_map = get_cell_map(grid)          # Mapping from reference to physical space
+cell_to_nodes = get_cell_node_ids(grid) # Node IDs for each cell
+node_coordinates = get_node_coordinates(grid) # Physical coordinates of nodes
+
+# ### Computing Cell-wise Node Coordinates
+#
+# There are two ways to get the coordinates of nodes for each cell:
+#
+# 1. Using standard Julia mapping:
+cell_to_node_coords = map(nodes -> node_coordinates[nodes], cell_to_nodes)
+
+# 2. Using Gridap's lazy evaluation system (more efficient for large meshes):
+cell_to_node_coords = lazy_map(Broadcasting(Reindex(node_coordinates)),cell_to_nodes)
+
+# ### Geometric Mappings
+#
+# The mapping from reference to physical space is defined by cell-wise linear combination of:
+# 1. Reference element shape functions (basis)
+# 2. Physical node coordinates (coefficients)
+
+cell_reffes = get_cell_reffe(grid)     # Get reference elements for each cell
+cell_basis = lazy_map(get_shapefuns,cell_reffes)  # Get basis functions
+cell_map = lazy_map(linear_combination,cell_to_node_coords,cell_basis)
+
+# ## 5. High-order Grids
+#
+# High-order geometric representations are essential for accurately modeling curved boundaries
+# and complex geometries. Let's explore this by creating a curved mesh:
+#
+# ### Example: Creating a Half-Cylinder
+#
+# First, we define a mapping that transforms our planar mesh into a half-cylinder:
+
+function F(x)
+  θ = x[1]*pi   # Map x-coordinate to angle [0,π]
+  z = x[2]      # Keep y-coordinate as height
+  VectorValue(cos(θ),sin(θ),z)  # Convert to cylindrical coordinates
+end
+
+# Apply the mapping to our node coordinates:
+new_node_coordinates = map(F,node_coordinates)
+
+# Create new cell mappings with the transformed coordinates:
+new_cell_to_node_coords = lazy_map(Broadcasting(Reindex(new_node_coordinates)),cell_to_nodes)
+new_cell_map = lazy_map(linear_combination,new_cell_to_node_coords,cell_basis)
+
+# Create a new grid with the transformed geometry:
+reffes, cell_types = compress_cell_data(cell_reffes)
+new_grid = UnstructuredGrid(new_node_coordinates,cell_to_nodes,reffes,cell_types)
+
+# Save for visualization:
+writevtk(new_grid,"half_cylinder_linear")
+
+#
+# If we visualize the result, we'll notice that despite applying a curved mapping,
+# our half-cylinder looks faceted. This is because we're still using linear elements
+# (straight edges) to approximate the curved geometry.
+#
+# ### Solution: High-order Elements
+#
+# To accurately represent curved geometries, we need high-order elements:
+
+# Create quadratic reference elements:
+order = 2  # Polynomial order
+new_reffes = [LagrangianRefFE(Float64,QUAD,order)]  # Quadratic quadrilateral elements
+new_cell_reffes = expand_cell_data(new_reffes,cell_types)
+
+# Create a finite element space to handle the high-order geometry:
+space = FESpace(model,new_cell_reffes)
+new_cell_to_nodes = get_cell_dof_ids(space)
+
+# Get the quadratic nodes in the reference element:
+cell_dofs = lazy_map(get_dof_basis,new_cell_reffes)
+cell_basis = lazy_map(get_shapefuns,new_cell_reffes)
+cell_to_ref_coordinates = lazy_map(get_nodes,cell_dofs)
+
+# Map the reference nodes to the physical space:
+cell_to_phys_coordinates = lazy_map(evaluate,cell_map,cell_to_ref_coordinates)
+
+# Create the high-order node coordinates:
+new_n_nodes = maximum(maximum,new_cell_to_nodes)
+new_node_coordinates = zeros(VectorValue{2,Float64},new_n_nodes)
+for (cell,nodes) in enumerate(new_cell_to_nodes)
+  for (i,node) in enumerate(nodes)
+    new_node_coordinates[node] = cell_to_phys_coordinates[cell][i]
+  end
+end
+
+# Apply our cylindrical mapping to the high-order nodes:
+new_node_coordinates = map(F,new_node_coordinates)
+
+# Create the high-order grid:
+new_grid = UnstructuredGrid(new_node_coordinates,new_cell_to_nodes,new_reffes,cell_types)
+writevtk(new_grid,"half_cylinder_quadratic")
+
+# The resulting mesh now accurately represents the curved geometry of the half-cylinder,
+# with quadratic elements properly capturing the curvature (despite paraview still showing 
+# a linear interpolation between the nodes).
+
+#
+# ## 6. Periodicity in Gridap
+#
+# Periodic boundary conditions are essential in many applications, such as:
+# - Modeling crystalline materials
+# - Simulating fluid flow in periodic domains
+# - Analyzing electromagnetic wave propagation
+#
+# Gridap handles periodicity through a clever approach:
+# 1. In the topology: Periodic vertices are "glued" together, creating a single topological entity
+# 2. In the geometry: The corresponding nodes maintain their distinct physical positions
+#
+# This separation between topology and geometry allows us to:
+# - Maintain the correct geometric representation
+# - Automatically enforce periodic boundary conditions
+# - Avoid mesh distortion at periodic boundaries
+#
+# ### Visualizing Periodicity
+#
+# Let's examine how periodicity affects the mesh structure through three examples:
+
+# 1. Standard non-periodic mesh:
+model = cartesian_model((0,1,0,1),(3,3))
+plot_node_numbering(model)
+# ![](../assets/geometry/nodes_nonperiodic.png)
+
+# 2. Mesh with periodicity in x-direction:
+model = cartesian_model((0,1,0,1),(3,3),isperiodic=(true,false))
+plot_node_numbering(model)
+# ![](../assets/geometry/nodes_halfperiodic.png)
+
+# 3. Mesh with periodicity in both directions:
+model = cartesian_model((0,1,0,1),(3,3),isperiodic=(true,true))
+plot_node_numbering(model)
+# ![](../assets/geometry/nodes_fullperiodic.png)
+
+# Notice how the vertex numbers (displayed at node positions) show the topological
+# connectivity, while the nodes remain at their physical positions.
+#
+# ### Creating a Möbius Strip
+#
+# We'll create it by:
+# 1. Starting with a rectangular mesh
+# 2. Making it periodic in one direction
+# 3. Adding a twist before connecting the ends
+#
+# #### Step 1: Create Base Mesh
+#
+# Start with a 3x3 cartesian mesh:
+
+nc = (3,3)  # Number of cells in each direction
+model = cartesian_model((0,1,0,1),nc)
+
+# Extract geometric and topological information:
+node_coords = get_node_coordinates(model)  # Physical positions
+cell_node_ids = get_cell_node_ids(model)  # Node connectivity
+cell_type = get_cell_type(model)          # Element type
+reffes = get_reffes(model)                # Reference elements
+
+# #### Step 2: Create Periodic Topology
+#
+# To create the Möbius strip, we need to:
+# 1. Identify vertices to be connected
+# 2. Reverse one edge to create the twist
+# 3. Ensure proper connectivity
+
+# Create initial vertex numbering:
+np = nc .+ 1  # Number of points in each direction
+mobius_ids = collect(LinearIndices(np))
+
+# Create the twist by reversing the last row:
+mobius_ids[end,:] = reverse(mobius_ids[1,:])
+
+# Map cell nodes to the new vertex numbering:
+cell_vertex_ids = map(nodes -> mobius_ids[nodes], cell_node_ids)
+
+# #### Step 3: Clean Up Vertex Numbering
+#
+# The new vertex IDs aren't contiguous (some numbers are duplicated due to periodicity).
+# We need to create a clean mapping:
+
+# Find unique vertices and create bidirectional mappings:
+vertex_to_node = unique(vcat(cell_vertex_ids...))
+node_to_vertex = find_inverse_index_map(vertex_to_node)
+
+# Renumber vertices to be contiguous:
+cell_vertex_ids = map(nodes -> node_to_vertex[nodes], cell_vertex_ids)
+
+# Convert to the Table format required by Gridap:
+cell_vertex_ids = Table(cell_vertex_ids)
+
+# Get coordinates for the vertices:
+vertex_coords = node_coords[vertex_to_node]
+polytopes = map(get_polytope,reffes)
+
+# #### Step 4: Create the Model
+#
+# Now we can create our Möbius strip model:
+
+# Create topology with periodic connections:
+topo = UnstructuredGridTopology(
+  vertex_coords, cell_vertex_ids, cell_type, polytopes
+)
+
+# Create grid (geometry remains unchanged):
+grid = UnstructuredGrid(
+  node_coords, cell_node_ids, reffes, cell_type
+)
+
+# Add basic face labels:
+labels = FaceLabeling(topo)
+
+# Combine into final model:
+mobius = UnstructuredDiscreteModel(grid,topo,labels)
+
+# Visualize the vertex numbering:
+plot_node_numbering(mobius)
+# ![](../assets/geometry/mobius.png)
diff --git a/src/poisson_amr.jl b/src/poisson_amr.jl
new file mode 100644
index 0000000..6d0deb1
--- /dev/null
+++ b/src/poisson_amr.jl
@@ -0,0 +1,198 @@
+#
+# In this tutorial, we will learn:
+#
+#    - How to use adaptive mesh refinement (AMR) with Gridap
+#    - How to set up a Poisson problem on an L-shaped domain
+#    - How to implement error estimation and marking strategies
+#    - How to visualize the AMR process and results
+#
+# ## Problem Overview
+#
+# We will solve the Poisson equation on an L-shaped domain using adaptive mesh refinement.
+# The L-shaped domain is known to introduce a singularity at its reentrant corner,
+# making it an excellent test case for AMR. The problem is:
+#
+# ```math
+# \begin{aligned}
+# -\Delta u &= f  &&\text{in } \Omega \\
+# u &= g &&\text{on } \partial\Omega
+# \end{aligned}
+# ```
+#
+# where Ω is the L-shaped domain, and we choose an exact solution with a singularity
+# to demonstrate the effectiveness of AMR.
+#
+# ## Error Estimation and AMR Process
+#
+# We use a residual-based a posteriori error estimator η. For each element K in the mesh,
+# the local error indicator ηK is computed as:
+#
+# ```math
+# \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
+# ```
+#
+# where:
+# - `h_K` is the diameter of element K
+# - `u_h` is the computed finite element solution
+# - The first term measures the element residual
+# - The second term measures the jump in the normal derivative across element boundaries
+#
+# The AMR process follows these steps in each iteration:
+#
+# 1. **Solve**: Compute the finite element solution u_h on the current mesh
+# 2. **Estimate**: Calculate error indicators ηK for each element
+# 3. **Mark**: Use Dörfler marking to select elements for refinement
+#    - Sort elements by error indicator
+#    - Mark elements containing a fixed fraction (here 80%) of total error
+# 4. **Refine**: Refine selected elements to obtain a new mesh. In this example, 
+#    we will be using the newest vertex bisection (NVB) method to keep the mesh 
+#    conforming (without any hanging nodes).
+#
+# This adaptive loop continues until either:
+# - A maximum number of iterations is reached
+# - The estimated error falls below a threshold
+# - The solution achieves desired accuracy
+#
+# The process automatically concentrates mesh refinement in regions of high error,
+# particularly around the reentrant corner where the solution has reduced regularity.
+# This results in better accuracy per degree of freedom compared to uniform refinement.
+#
+# ## Required Packages
+
+using Gridap, Gridap.Geometry, Gridap.Adaptivity
+using DataStructures
+
+# ## Problem Setup
+#
+# We define an exact solution that contains a singularity at the corner (0.5, 0.5)
+# of the L-shaped domain. This singularity will demonstrate how AMR automatically
+# refines the mesh in regions of high error.
+
+ϵ = 1e-2
+r(x) = ((x[1]-0.5)^2 + (x[2]-0.5)^2)^(1/2)
+u_exact(x) = 1.0 / (ϵ + r(x))
+
+# Create an L-shaped domain by removing a quadrant from a unit square.
+# The domain is [0,1]² \ [0.5,1]×[0.5,1]
+function LShapedModel(n)
+  model = CartesianDiscreteModel((0,1,0,1),(n,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))
+end
+
+# Define the L2 norm for error estimation.
+# These will be used to compute both local and global error measures.
+l2_norm(he,xh,dΩ) = ∫(he*(xh*xh))*dΩ
+l2_norm(xh,dΩ) = ∫(xh*xh)*dΩ
+
+# ## AMR Step Function
+#
+# The `amr_step` function performs a single step of the adaptive mesh refinement process:
+# 1. Solves the Poisson problem on the current mesh
+# 2. Estimates the error using residual-based error indicators
+# 3. Marks cells for refinement using Dörfler marking
+# 4. Refines the mesh using newest vertex bisection (NVB)
+
+function amr_step(model,u_exact;order=1)
+  "Create FE spaces with Dirichlet boundary conditions on all boundaries"
+  reffe = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(model,reffe;dirichlet_tags=["boundary"])
+  U = TrialFESpace(V,u_exact)
+  
+  "Setup integration measures"
+  Ω = Triangulation(model)
+  Γ = Boundary(model)
+  Λ = Skeleton(model)
+  
+  dΩ = Measure(Ω,4*order)
+  dΓ = Measure(Γ,2*order)
+  dΛ = Measure(Λ,2*order)
+  
+  "Compute cell sizes for error estimation"
+  hK = CellField(sqrt.(collect(get_array(∫(1)dΩ))),Ω)
+
+  "Get normal vectors for boundary and interface terms"
+  nΓ = get_normal_vector(Γ)
+  nΛ = get_normal_vector(Λ)
+
+  "Define the weak form"
+  ∇u(x)  = ∇(u_exact)(x)
+  f(x)   = -Δ(u_exact)(x)
+  a(u,v) = ∫(∇(u)⋅∇(v))dΩ
+  l(v)   = ∫(f*v)dΩ
+  
+  "Define the residual error estimator
+  It includes volume residual, boundary jump, and interface jump terms"
+  ηh(u)  = l2_norm(hK*(f + Δ(u)),dΩ) +           # Volume residual
+           l2_norm(hK*(∇(u) - ∇u)⋅nΓ,dΓ) +       # Boundary jump
+           l2_norm(jump(hK*∇(u)⋅nΛ),dΛ)          # Interface jump
+  
+  "Solve the FE problem"
+  op = AffineFEOperator(a,l,U,V)
+  uh = solve(op)
+  
+  "Compute error indicators"
+  η = estimate(ηh,uh)
+  
+  "Mark cells for refinement using Dörfler marking
+  This strategy marks cells containing a fixed fraction (0.9) of the total error"
+  m = DorflerMarking(0.9)
+  I = Adaptivity.mark(m,η)
+  
+  "Refine the mesh using newest vertex bisection"
+  method = Adaptivity.NVBRefinement(model)
+  amodel = refine(method,model;cells_to_refine=I)
+  fmodel = Adaptivity.get_model(amodel)
+
+  "Compute the global error for convergence testing"
+  error = sum(l2_norm(uh - u_exact,dΩ))
+  return fmodel, uh, η, I, error
+end
+
+# ## Main AMR Loop
+#
+# We perform multiple AMR steps, refining the mesh iteratively and solving
+# the problem on each refined mesh. This demonstrates how the error decreases
+# as the mesh is adaptively refined in regions of high error.
+
+nsteps = 5
+order = 1
+model = LShapedModel(10)
+
+last_error = Inf
+for i in 1:nsteps
+  fmodel, uh, η, I, error = amr_step(model,u_exact;order)
+  
+  is_refined = map(i -> ifelse(i ∈ I, 1, -1), 1:num_cells(model))
+  
+  Ω = Triangulation(model)
+  writevtk(
+    Ω,"model_$(i-1)",append=false,
+    cellfields = [
+      "uh" => uh,                    # Computed solution
+      "η" => CellField(η,Ω),        # Error indicators
+      "is_refined" => CellField(is_refined,Ω),  # Refinement markers
+      "u_exact" => CellField(u_exact,Ω),       # Exact solution
+    ],
+  )
+  
+  println("Error: $error, Error η: $(sum(η))")
+  last_error = error
+  model = fmodel
+end
+
+# The final mesh gives the following result:
+# ![](../assets/amr/result.png)
+#
+# ## Conclusion
+#
+# In this tutorial, we have demonstrated how to:
+# 1. Implement adaptive mesh refinement for the Poisson equation
+# 2. Use residual-based error estimation to identify regions needing refinement
+# 3. Apply Dörfler marking to select cells for refinement
+# 4. Visualize the AMR process and solution convergence
+#
+# The results show how AMR automatically refines the mesh near the singularity,
+# leading to more efficient and accurate solutions compared to uniform refinement.
diff --git a/src/unfitted_poisson.jl b/src/poisson_unfitted.jl
similarity index 100%
rename from src/unfitted_poisson.jl
rename to src/poisson_unfitted.jl