diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0d3a5f4..7c3c8d8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # CHANGES
 
+## v1.6.0
+  - added Curl2D and Curl3D evaluation for Hdiv
+  
+
 ## v1.5.0 November 21, 2025
   - added type parameter to `PointEvaluator` so `lazy_interpolate!` is now accessible for automatic differentiation
   - added new `compute_lazy_interpolation_jacobian` function for computing Jacobians of interpolations between finite element spaces
diff --git a/Project.toml b/Project.toml
index 2698e2e..aea5af2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ExtendableFEMBase"
 uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
-version = "1.5.1"
+version = "1.6.0"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
 
 [deps]
diff --git a/src/feevaluator_hdiv.jl b/src/feevaluator_hdiv.jl
index 8882522..6419842 100644
--- a/src/feevaluator_hdiv.jl
+++ b/src/feevaluator_hdiv.jl
@@ -106,3 +106,49 @@ function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Grad
     end
     return nothing
 end
+
+
+# CURL2D HDIV
+function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl2D, <:AbstractHdivFiniteElement})
+    L2GAinv = _update_trafo!(FEBE)
+    L2GM = _update_piola!(FEBE)
+    subset = _update_subset!(FEBE)
+    coefficients = _update_coefficients!(FEBE)
+    cvals = FEBE.cvals
+    offsets2 = FEBE.offsets2
+    refbasisderivvals = FEBE.refbasisderivvals
+    fill!(cvals, 0)
+    det = FEBE.L2G.det # 1 alloc
+    for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2)
+        for j in 1:size(L2GM, 2), m in 1:size(L2GAinv, 2)
+            cvals[1, dof_i, i] -= L2GAinv[2, m] * L2GM[1, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[1, dof_i] / det
+            cvals[1, dof_i, i] += L2GAinv[1, m] * L2GM[2, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[2, dof_i] / det
+        end
+    end
+    return nothing
+end
+
+
+# CURL3D HDIV
+function update_basis!(FEBE::SingleFEEvaluator{<:Real, <:Real, <:Integer, <:Curl3D, <:AbstractHdivFiniteElement})
+    L2GAinv = _update_trafo!(FEBE)
+    L2GM = _update_piola!(FEBE)
+    subset = _update_subset!(FEBE)
+    coefficients = _update_coefficients!(FEBE)
+    cvals = FEBE.cvals
+    offsets2 = FEBE.offsets2
+    refbasisderivvals = FEBE.refbasisderivvals
+    fill!(cvals, 0)
+    det = FEBE.L2G.det # 1 alloc
+    for i in 1:size(cvals, 3), dof_i in 1:size(cvals, 2)
+        for j in 1:size(L2GM, 2), m in 1:size(L2GAinv, 2)
+            cvals[1, dof_i, i] -= L2GAinv[3, m] * L2GM[2, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[2, dof_i] / det
+            cvals[1, dof_i, i] += L2GAinv[2, m] * L2GM[3, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[3, dof_i] / det
+            cvals[2, dof_i, i] -= L2GAinv[1, m] * L2GM[3, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[3, dof_i] / det
+            cvals[2, dof_i, i] += L2GAinv[3, m] * L2GM[1, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[1, dof_i] / det
+            cvals[3, dof_i, i] -= L2GAinv[2, m] * L2GM[1, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[1, dof_i] / det
+            cvals[3, dof_i, i] += L2GAinv[1, m] * L2GM[2, j] * refbasisderivvals[subset[dof_i] + offsets2[j], m, i] * coefficients[2, dof_i] / det
+        end
+    end
+    return nothing
+end
diff --git a/test/test_operators.jl b/test/test_operators.jl
index bafe31a..a943e8e 100644
--- a/test/test_operators.jl
+++ b/test/test_operators.jl
@@ -6,8 +6,12 @@ function run_operator_tests()
         println("============================")
         error = test_derivatives2D()
         @test error < 1.0e-14
+        error = test_derivatives2D_hdiv()
+        @test error < 1.0e-14
         error = test_derivatives3D()
         @test error < 1.0e-14
+        error = test_derivatives3D_hdiv()
+        @test error < 1.0e-14
         test_reconstructions()
     end
 end
@@ -79,8 +83,7 @@ function test_derivatives2D()
     xgrid = grid_triangle([-1.0 0.0; 1.0 0.0; 0.0 1.0]')
 
     ## define P2-Courant finite element space
-    FEType = H1P2{2, 2}
-    FES = FESpace{FEType}(xgrid)
+    FES = FESpace{H1P2{2, 2}}(xgrid)
     show(devnull, FES)
 
     ## get midpoint quadrature rule for constants
@@ -140,6 +143,52 @@ function test_derivatives2D()
 end
 
 
+function test_derivatives2D_hdiv()
+    ## define test function and expected operator evals
+    function testf(result, qpinfo)
+        x = qpinfo.x
+        result[1] = x[1]^2
+        return result[2] = 3 * x[2]^2 + x[1] * x[2]
+    end
+
+    ## expected values of operators in cell midpoint
+    expected_curl2 = [1 / 3]
+
+    ## define grid = a single non-refenrece triangle
+    xgrid = grid_triangle([-1.0 0.0; 1.0 0.0; 0.0 1.0]')
+
+    ## define P2-Courant finite element space
+    FES = FESpace{HDIVBDM2{2}}(xgrid)
+    show(devnull, FES)
+
+    ## get midpoint quadrature rule for constants
+    qf = QuadratureRule{Float64, Triangle2D}(0)
+
+    ## define FE basis Evaluator for Hessian
+    FEBE_curl2 = FEEvaluator(FES, Curl2D, qf)
+
+    ## update on cell 1
+    update_basis!(FEBE_curl2, 1)
+
+    ## interpolate quadratic testfunction
+    Iu = FEVector(FES)
+    interpolate!(Iu[1], testf)
+
+    ## check if operator evals have the correct length
+    @assert size(FEBE_curl2.cvals, 1) == length(expected_curl2)
+
+    # compute derivatives
+    curl2 = zeros(Float64, 1)
+    eval_febe!(curl2, FEBE_curl2, Iu.entries[FES[CellDofs][:, 1]], 1)
+
+    ## compute errors to expected values
+    error_curl2 = sqrt(sum((curl2 - expected_curl2) .^ 2))
+    println("EG = Triangle2D | operator = Curl2 | error = $error_curl2")
+
+    return maximum([error_curl2])
+end
+
+
 function test_derivatives3D()
     ## define test function and expected operator evals
     function testf(result, qpinfo)
@@ -220,3 +269,53 @@ function test_derivatives3D()
 
     return maximum([error_curl3, error_L, error_H, error_symH, error_symH2])
 end
+
+
+function test_derivatives3D_hdiv()
+    ## define test function and expected operator evals
+    function testf(result, qpinfo)
+        x = qpinfo.x
+        result[1] = x[2]
+        result[2] = 3 * x[3]
+        return result[3] = x[2]
+    end
+
+    ## expected values of operators in cell midpoint
+    expected_curl3 = [1 - 3, 0, -1]
+
+    ## define grid = a single non-refenrece triangle
+    xgrid = reference_domain(Tetrahedron3D)
+    xgrid[Coordinates][:, 2] = [2, 0, 0]
+
+    ## define P2-Courant finite element space
+    FES = FESpace{HDIVRT1{3}}(xgrid)
+    show(devnull, FES)
+
+    ## get midpoint quadrature rule for constants
+    qf = QuadratureRule{Float64, Tetrahedron3D}(0)
+
+    ## define FE basis Evaluator for Hessian
+    FEBE_curl3 = FEEvaluator(FES, Curl3D, qf)
+
+    ## update on cell 1
+    update_basis!(FEBE_curl3, 1)
+
+    ## interpolate quadratic testfunction
+    Iu = FEVector(FES)
+    interpolate!(Iu[1], testf)
+
+    ## check if operator evals have the correct length
+    @assert size(FEBE_curl3.cvals, 1) == length(expected_curl3)
+
+    ## eval 2nd order derivatives at only quadrature point 1
+    ## since function is quadratic this should be constant
+    curl3 = zeros(Float64, 3)
+    eval_febe!(curl3, FEBE_curl3, Iu.entries[FES[CellDofs][:, 1]], 1)
+    @info curl3
+
+    ## compute errors to expected values
+    error_curl3 = sqrt(sum((curl3 - expected_curl3) .^ 2))
+    println("EG = Tetrahedron3D | operator = Curl3 | error = $error_curl3")
+
+    return maximum([error_curl3])
+end