diff --git a/NEWS.md b/NEWS.md
index e9287dbdc..524a03951 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,12 @@ 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]
+
+### Added
+
+- Added AMR-related methods `mark` and `estimate` to `Adaptivity` module. Implemented Dorfler marking strategy. Since PR[#1063](https://github.com/gridap/Gridap.jl/pull/1063).
+
 ## [0.18.8] - 2024-12-2
 
 ### Added
diff --git a/docs/src/Adaptivity.md b/docs/src/Adaptivity.md
index 377cd5d23..5a698528a 100644
--- a/docs/src/Adaptivity.md
+++ b/docs/src/Adaptivity.md
@@ -80,6 +80,18 @@ The API is given by the following methods:
   MacroReferenceFE
 ```
 
+## Adaptive Mesh Refinement
+
+One of the main uses of mesh refinement is Adaptive Mesh Refinement, where the mesh is refined only in regions of interest.
+
+The typical AMR workflow is the so-called `solve-estimate-mark-refine` loop. Since estimators will generally be problem-dependent, we only aim to provide some generic tools that can be combined by the user:
+
+```@docs
+  DorflerMarking
+  mark
+  estimate
+```
+
 ## Notes for users
 
 Most of the tools provided by this module are showcased in the tests of the module itself, as well as the following tutorial (coming soon).
diff --git a/src/Adaptivity/AdaptiveMeshRefinement.jl b/src/Adaptivity/AdaptiveMeshRefinement.jl
new file mode 100644
index 000000000..40f9eea42
--- /dev/null
+++ b/src/Adaptivity/AdaptiveMeshRefinement.jl
@@ -0,0 +1,188 @@
+
+# Marking strategies
+
+"""
+    struct DorflerMarking
+      θ :: Float64
+      ν :: Float64
+      strategy :: Symbol
+    end
+
+    DorflerMarking(θ::Float64; ν::Float64 = 0.5, strategy::Symbol = :quickmark)
+
+Implements the Dorfler marking strategy. Given a vector `η` of real positive numbers, 
+the marking strategy find a subset of indices `I` such that 
+
+  sum(η[I]) > θ * sum(η)
+
+where `0 < θ < 1` is a threshold parameter. 
+
+For more details, see the following reference: 
+
+"Dörfler marking with minimal cardinality is a linear complexity problem", Pfeiler et al. (2020)
+
+The marking algorithm is controlled by the `strategy` parameter, which can take 
+the following values: 
+
+- `:sort`: Optimal cardinality, O(N log N) complexity. See Algorithm 2 in the reference.
+- `:binsort`: Quasi-optimal cardinality, O(N) complexity. See Algorithm 7 in the reference.
+- `:quickmark`: Optimal cardinality, O(N) complexity.  See Algorithm 10 in the reference.
+
+# Arguments
+
+- `θ::Float64`: The threshold parameter. Between 0 and 1.
+- `ν::Float64`: Extra parameter for `:binsort`. Default is 0.5.
+- `strategy::Symbol`: The marking strategy. Default is `:quickmark`.
+
+# Usage
+
+```julia
+η = abs.(randn(1000))
+m = DorflerMarking(0.5)
+I = mark(m,η)
+```
+
+"""
+struct DorflerMarking
+  θ :: Float64
+  ν :: Float64
+  strategy :: Symbol
+  function DorflerMarking(
+    θ::Float64;
+    ν::Float64 = 0.5,
+    strategy::Symbol = :quickmark
+  )
+    @assert 0 < θ < 1
+    @assert strategy ∈ (:sort,:binsort,:quickmark) "Strategy not recognized. Available values are (:sort)"
+    new(θ,ν,strategy)
+  end
+end
+
+"""
+    mark(m::DorflerMarking, η::Vector{<:Real}) -> Vector{Int}
+
+Given a vector `η` of real positive numbers, returns a subset of indices `I` such that 
+satisfying the Dorfler marking condition.
+"""
+mark(m::DorflerMarking, η::Vector{<:Real}) = mark(Val(m.strategy), m, η)
+
+function mark(::Val{:sort}, m::DorflerMarking, η::Vector{<:Real})
+  target = m.θ * sum(η)
+  perm = sortperm(η, rev=true, alg=QuickSort)
+  s = zero(eltype(η))
+  k = 0
+  while s < target
+    k += 1
+    s += η[perm[k]]
+  end
+  return perm[1:k]
+end
+
+function mark(::Val{:binsort}, m::DorflerMarking, η::Vector{<:Real})
+  target = m.θ * sum(η)
+  M = maximum(η)
+  N = length(η)
+  
+  # Find minimal K such that
+  #   νᴷ⁺¹ M ≤ (1 - θ) / (θ * N) * target
+  K = 0
+  a = M*m.ν
+  b = (1 - m.θ) / (m.θ * N * M) * target
+  while a > b
+    K += 1
+    a *= m.ν
+  end
+
+  # Sort into bins Bk such that ηi ∈ Bk if
+  #    νᴷ⁺¹ M ≤ ηi < νᴷ M
+  bins = zeros(Int,N)
+  sums = zeros(Float64,K+1)
+  lbs = zeros(Float64,K+1)
+  lbs[1] = M * m.ν
+  for i in 2:K
+    lbs[i] = lbs[i-1] * m.ν
+  end
+  lbs[end] = 0.0
+
+  for (i,ηi) in enumerate(η)
+    k = 1
+    while ηi < lbs[k]
+      k += 1
+    end
+    bins[i] = k
+    sums[k] += ηi
+  end
+
+  # Find minimal set of bins that gets over target
+  k = 0
+  s = 0.0
+  while s < target
+    k += 1
+    s += sums[k]
+  end
+  
+  return findall(i -> i <= k, bins)
+end
+
+function mark(::Val{:quickmark}, m::DorflerMarking, η::Vector{<:Real})
+  function quickmark!(η, perm, l, u, target) :: Int
+    m = (u - l) ÷ 2
+    sort!(view(perm,l:u), by=i->η[i], rev=true, alg=PartialQuickSort(m))
+
+    p = l + m
+    t = η[perm[p]]
+    σ = sum(η[perm[l:p-1]])
+    
+    (σ >= target) && return quickmark!(η, perm, l, p, target)
+    (σ + t >= target) && return p
+    return quickmark!(η, perm, p + 1, u, target - σ - t)
+  end
+
+  N = length(η)
+  l, u = 1, N
+  perm = collect(1:N)
+  target = m.θ * sum(η)
+  m = quickmark!(η, perm, l, u, target)
+  return perm[1:m]
+end
+
+# Estimators
+
+"""
+    estimate(f::Function, uh::Function) -> Vector{Float64}
+
+Given a functional `f` and a function `uh`, such that `f(uh)` produces a 
+scalar-valued `DomainContribution`, collects the estimator values for 
+each cell in the background model.
+"""
+function estimate(f::Function, uh)
+  collect_estimator(f(uh))
+end
+
+function collect_estimator(c::DomainContribution)
+  trians = get_domains(c)
+  bgmodel = get_background_model(first(trians))
+  msg = "Estimator not implemented for mixed background models"
+  @notimplementedif !all([bgmodel == get_background_model(trian) for trian in trians]) msg
+
+  Dc = num_cell_dims(bgmodel)
+  η = zeros(Float64,num_cells(bgmodel))
+  for trian in trians
+    glue = get_glue(trian,Val(Dc))
+    collect_estimator!(η,glue,get_contribution(c,trian))
+  end
+
+  return η
+end
+
+function collect_estimator!(η, glue::FaceToFaceGlue, c)
+  cache = array_cache(c)
+  for (face,bgcell) in enumerate(glue.tface_to_mface)
+    η[bgcell] += getindex!(cache,c,face)
+  end
+end
+
+function collect_estimator!(η, glue::SkeletonPair, c)
+  collect_estimator!(η,glue.plus,c)
+  collect_estimator!(η,glue.minus,c)
+end
diff --git a/src/Adaptivity/Adaptivity.jl b/src/Adaptivity/Adaptivity.jl
index d80eb0c3e..32a5d1594 100644
--- a/src/Adaptivity/Adaptivity.jl
+++ b/src/Adaptivity/Adaptivity.jl
@@ -40,6 +40,8 @@ export AdaptedTriangulation
 export Triangulation, is_change_possible, best_target, get_adapted_model
 export change_domain, move_contributions
 
+export DorflerMarking, mark, estimate
+
 include("RefinementRules.jl")
 include("FineToCoarseFields.jl")
 include("OldToNewFields.jl")
@@ -51,5 +53,6 @@ include("MacroFEs.jl")
 include("CompositeQuadratures.jl")
 include("EdgeBasedRefinement.jl")
 include("SimplexifyRefinement.jl")
+include("AdaptiveMeshRefinement.jl")
 
 end # module
diff --git a/src/Adaptivity/EdgeBasedRefinement.jl b/src/Adaptivity/EdgeBasedRefinement.jl
index 724f74fd4..fc50d45ea 100644
--- a/src/Adaptivity/EdgeBasedRefinement.jl
+++ b/src/Adaptivity/EdgeBasedRefinement.jl
@@ -78,7 +78,6 @@ function refine_edge_based_topology(
   c2n_map_new = get_refined_cell_to_vertex_map(topo,rrules,faces_list)
   polys_new, cell_type_new = _get_cell_polytopes(rrules)
   orientation = NonOriented()
-
   return UnstructuredGridTopology(coords_new,c2n_map_new,cell_type_new,polys_new,orientation)
 end
 
@@ -109,17 +108,28 @@ The new vertices are ordered by parent dimension and face id (in that order).
 function get_new_coordinates_from_faces(p::Union{Polytope{D},GridTopology{D}},faces_list::Tuple) where {D}
   @check length(faces_list) == D+1
 
-  nN_new     = sum(x->length(x),faces_list)
+  nN_new     = sum(length,faces_list)
   coords_old = get_vertex_coordinates(p)
   coords_new = Vector{eltype(coords_old)}(undef,nN_new)
 
   n = 1
-  for (d,dfaces) in enumerate(faces_list)
-    if length(dfaces) > 0
-      nf = length(dfaces)
-      d2n_map = get_faces(p,d-1,0)
-      coords_new[n:n+nf-1] .= map(f -> sum(coords_old[d2n_map[f]])/length(d2n_map[f]), dfaces)
-      n += nf
+  # Nodes
+  if !isempty(faces_list[1])
+    for node in faces_list[1]
+      coords_new[n] = coords_old[node]
+      n += 1
+    end
+  end
+  # Faces (d > 0)
+  for (d,dfaces) in enumerate(faces_list[2:end])
+    if !isempty(dfaces)
+      d2n_map = get_faces(p,d,0)
+      cache = array_cache(d2n_map)
+      for face in dfaces
+        face_nodes = getindex!(cache,d2n_map,face)
+        coords_new[n] = sum(coords_old[face_nodes])/length(face_nodes)
+        n += 1
+      end
     end
   end
 
@@ -255,12 +265,13 @@ function setup_edge_based_rrules(method::NVBRefinement, topo::UnstructuredGridTo
   setup_edge_based_rrules(method, topo, collect(1:num_faces(topo,Dc)))
 end
 
-function setup_edge_based_rrules(method::NVBRefinement, topo::UnstructuredGridTopology{Dc},cells_to_refine::AbstractArray{<:Integer}) where Dc
+function setup_edge_based_rrules(
+  method::NVBRefinement, topo::UnstructuredGridTopology{Dc}, cells_to_refine::AbstractArray{<:Integer}
+) where Dc
   nE = num_faces(topo,1)
-  c2e_map       = get_faces(topo,Dc,1)
-  c2e_map_cache = array_cache(c2e_map)
-  e2c_map       = get_faces(topo,1,Dc)
-  polys       = topo.polytopes
+  c2e_map = get_faces(topo,Dc,1)
+  e2c_map = get_faces(topo,1,Dc)
+  polys   = topo.polytopes
   cell_types  = topo.cell_type
   cell_color  = copy(cell_types) # WHITE
   # Hardcoded for TRI
@@ -274,43 +285,39 @@ function setup_edge_based_rrules(method::NVBRefinement, topo::UnstructuredGridTo
   c_to_longest_edge_lid = method.cell_to_longest_edge_lid
   is_refined = falses(nE)
   # Loop over cells and mark edges to refine i.e. is_refined
-  # The reason to not loop directly on c is that we need to change c within
-  # a given iteration of the for loop
-  for i in 1:length(cells_to_refine)
-    c = cells_to_refine[i]
-    e_longest = c_to_longest_edge_gid[c]
+  for i in eachindex(cells_to_refine)
     # Has to terminate because an edge is marked each iteration or we skip an
     # iteration due to a boundary cell
-    while !is_refined[e_longest]
-      is_refined[e_longest] = true
-      c_nbor_lid = findfirst(c′ -> c′ != c, e2c_map[e_longest])
-      if isnothing(c_nbor_lid) # We've reach the boundary
+    c = cells_to_refine[i]
+    e = c_to_longest_edge_gid[c]
+    while !is_refined[e]
+      is_refined[e] = true
+      e_cells = view(e2c_map,e)
+      if length(e_cells) == 1 # We've reach the boundary
         continue
       else
-        # Get the longest edge of the neighbor
-        c_nbor_gid = e2c_map[e_longest][c_nbor_lid]
-        e_longest = c_to_longest_edge_gid[c_nbor_gid]
-        # Set the current cell gid to that of the neighbor
-        c = c_nbor_gid
+        # Propagate to neighboring cell
+        c = ifelse(e_cells[1] == c, e_cells[2], e_cells[1])
+        e = c_to_longest_edge_gid[c]
       end
     end
   end
   # Loop over cells and refine based on marked edges
   for c in 1:length(c2e_map)
-    c_edges = getindex!(c2e_map_cache, c2e_map, c)
+    c_edges = view(c2e_map,c)
     refined_edge_lids = findall(is_refined[c_edges])
-    # GREEN refinement because only one edge should be bisected
     if length(refined_edge_lids) == 1
+      # GREEN refinement because only one edge should be bisected
       ref_edge = refined_edge_lids[1]
       cell_color[c] = GREEN + Int8(ref_edge-1)
-      # BLUE refinement: two bisected
     elseif length(refined_edge_lids) == 2
+      # BLUE refinement: two bisected
       long_ref_edge_lid = c_to_longest_edge_lid[c]
       short_ref_edge_lid = setdiff(refined_edge_lids, long_ref_edge_lid)[1]
       blue_idx = BLUE_dict[(long_ref_edge_lid, short_ref_edge_lid)]
       cell_color[c] = BLUE + Int8(blue_idx - 1)
-      # DOUBLE BLUE refinement: three bisected edges (somewhat rare)
     elseif length(refined_edge_lids) == 3
+      # DOUBLE BLUE refinement: three bisected edges (somewhat rare)
       long_ref_edge_lid = c_to_longest_edge_lid[c]
       cell_color[c] = BLUE_DOUBLE + Int(long_ref_edge_lid - 1)
     end
diff --git a/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl
new file mode 100644
index 000000000..6f0915d0e
--- /dev/null
+++ b/test/AdaptivityTests/AdaptiveMeshRefinementTests.jl
@@ -0,0 +1,112 @@
+module AdaptiveMeshRefinementTests
+
+using Test
+using Gridap, Gridap.Geometry, Gridap.Adaptivity
+using DataStructures
+
+# Marking tests
+
+function test_dorfler_marking()
+  for strategy in (:sort,:binsort,:quickmark)
+    for n in (1000,10000)
+      for θ in (0.3,0.5)
+        η = abs.(randn(n))
+        m = DorflerMarking(θ;strategy)
+        I = Adaptivity.mark(m,η)
+        @test sum(η[I]) > θ * sum(η)
+        println("Strategy: $strategy, θ: $θ, n: $n")
+        println("  > N marked = $(length(I)), val marked = $(sum(η[I]) / sum(η))")
+      end
+    end
+  end
+end
+
+# AMR tests
+
+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
+
+l2_norm(he,xh,dΩ) = ∫(he*(xh*xh))*dΩ
+l2_norm(xh,dΩ) = ∫(xh*xh)*dΩ
+
+function amr_step(model,u_exact;order=1)
+  reffe = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(model,reffe;dirichlet_tags=["boundary"])
+  U = TrialFESpace(V,u_exact)
+  
+  Ω = Triangulation(model)
+  Γ = Boundary(model)
+  Λ = Skeleton(model)
+  
+  dΩ = Measure(Ω,4*order)
+  dΓ = Measure(Γ,2*order)
+  dΛ = Measure(Λ,2*order)
+  
+  hK = CellField(sqrt.(collect(get_array(∫(1)dΩ))),Ω)
+
+  nΓ = get_normal_vector(Γ)
+  nΛ = get_normal_vector(Λ)
+
+  ∇u(x)  = ∇(u_exact)(x)
+  f(x)   = -Δ(u_exact)(x)
+  a(u,v) = ∫(∇(u)⋅∇(v))dΩ
+  l(v)   = ∫(f*v)dΩ
+  ηh(u)  = l2_norm(hK*(f + Δ(u)),dΩ) + l2_norm(hK*(∇(u) - ∇u)⋅nΓ,dΓ) + l2_norm(jump(hK*∇(u)⋅nΛ),dΛ)
+  
+  op = AffineFEOperator(a,l,U,V)
+  uh = solve(op)
+  η = estimate(ηh,uh)
+  
+  m = DorflerMarking(0.8)
+  I = Adaptivity.mark(m,η)
+  
+  method = Adaptivity.NVBRefinement(model)
+  fmodel = Adaptivity.get_model(refine(method,model;cells_to_refine=I))
+
+  error = sum(l2_norm(uh - u_exact,dΩ))
+  return fmodel, uh, η, I, error
+end
+
+function test_amr(nsteps,order)
+  model = LShapedModel(10)
+
+  ϵ = 1e-2
+  r(x) = ((x[1]-0.5)^2 + (x[2]-0.5)^2)^(1/2)
+  u_exact(x) = 1.0 / (ϵ + r(x))
+
+  vtk = false
+  last_error = Inf
+  for i in 1:nsteps
+    fmodel, uh, η, I, error = amr_step(model,u_exact;order)
+    if vtk
+      is_refined = map(i -> ifelse(i ∈ I, 1, -1), 1:num_cells(model))
+      Ω = Triangulation(model)
+      writevtk(
+        Ω,"tmp/model_$(i-1)",append=false,
+        cellfields = [
+          "uh" => uh,
+          "η" => CellField(η,Ω),
+          "is_refined" => CellField(is_refined,Ω),
+          "u_exact" => CellField(u_exact,Ω),
+        ],
+      )
+    end
+
+    println("Error: $error, Error η: $(sum(η))")
+    @test (i < 3) || (error < last_error)
+    last_error = error
+    model = fmodel
+  end
+end
+
+############################################################################################
+
+@testset "Dorfler marking" test_dorfler_marking()
+@testset "AMR - Poisson" test_amr(20,2)
+
+end # module
\ No newline at end of file
diff --git a/test/AdaptivityTests/runtests.jl b/test/AdaptivityTests/runtests.jl
index 87437afc3..2d2a569da 100644
--- a/test/AdaptivityTests/runtests.jl
+++ b/test/AdaptivityTests/runtests.jl
@@ -23,4 +23,8 @@ end
   include("MacroFEStokesTests.jl")
 end
 
+@testset "AMR" begin
+  include("AdaptiveMeshRefinementTests.jl")
+end
+
 end # module
\ No newline at end of file