diff --git a/Project.toml b/Project.toml
index 221cab7..d29fb2a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ReliabilityOptimization"
 uuid = "55eddd50-7f45-4398-96cf-6a37e2b16f80"
 authors = ["Lucas Pereira and Mohamed Tarek"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -27,14 +27,18 @@ FiniteDifferences = "0.12"
 ImplicitDifferentiation = "0.2"
 JuliaFormatter = "1"
 NonconvexIpopt = "0.4"
+NonconvexTOBS = "0.2"
 Reexport = "1"
 StaticArraysCore = "1"
+TopOpt = "0.8"
 UnPack = "1"
 Zygote = "0.6"
 julia = "1"
 
 [extras]
+NonconvexTOBS = "6c0b5230-d4c9-466e-bfd4-b31e6272ab65"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TopOpt = "53a1e1a5-51bb-58a9-8a02-02056cc81109"
 
 [targets]
-test = ["Test"]
+test = ["NonconvexTOBS", "Test", "TopOpt"]
diff --git a/src/ReliabilityOptimization.jl b/src/ReliabilityOptimization.jl
index 4188f63..ce249df 100644
--- a/src/ReliabilityOptimization.jl
+++ b/src/ReliabilityOptimization.jl
@@ -101,14 +101,12 @@ function (f::RandomFunction)(x)
     return MvNormal(muf, covf)
 end
 
-# necessary type piracy FiniteDifferences._estimate_magnitudes uses this constructor which Zygote struggles to differentiate on its own
-function ChainRulesCore.rrule(
-    ::typeof(StaticArraysCore.SVector{3}),
-    x1::T,
-    x2::T,
-    x3::T,
-) where {T}
-    StaticArraysCore.SVector{3}(x1, x2, x3), Δ -> (NoTangent(), Δ[1], Δ[2], Δ[3])
+function ChainRulesCore.ProjectTo(x::StaticArraysCore.SVector{N, Vector{Float64}}) where {N}
+    return ChainRulesCore.ProjectTo{SArray}(;
+        element = ChainRulesCore.ProjectTo(zero(first(x))),
+        axes = axes(x),
+        size = Size(x),
+    )
 end
 
 end
diff --git a/test/example.jl b/test/example.jl
index 19ed83f..4e10cd0 100644
--- a/test/example.jl
+++ b/test/example.jl
@@ -26,7 +26,7 @@ function uncertainComp(x, logEs)
     Es = exp.(logEs)
     # interpolation of properties between materials
     interp = MaterialInterpolation(Es, penalty)
-    MultiMaterialVariables(x, nmats) |> interp |> filter |> comp
+    MultiMaterialVariables(x, nmats) |> tounit |> filter |> interp |> comp
     # return sum(x) + sum(Es)
 end
 # wrap original function in RandomFunction struct
@@ -37,19 +37,29 @@ x0 = fill(M, ncells * (length(logEs) - 1))
 # (Returns propability distribution of the objective for current point)
 d = rf(x0)
 # mass constraint
+constr_penalty = TopOpt.PowerPenalty(1.0)
+constr_interp = MaterialInterpolation(densities, constr_penalty)
 constr = x -> begin
-    ρs = PseudoDensities(MultiMaterialVariables(x, nmats))
-    return sum(element_densities(ρs, densities)) / ncells - 0.3 # unit element volume
+    ρs = constr_interp(MultiMaterialVariables(x, nmats))
+    return sum(ρs.x) / ncells - 0.3 # elements have unit volumes
 end
 function obj(x) # objective for TO problem
     dist = rf(x)
     mean(dist)[1] + 2 * sqrt(cov(dist)[1, 1])
 end
 obj(x0)
-Zygote.gradient(obj, x0)
-FiniteDifferences.grad(central_fdm(5, 1), obj, x0)[1]
+g1 = Zygote.gradient(obj, x0)[1]
+g2 = FiniteDifferences.grad(central_fdm(5, 1), obj, x0)[1]
+norm(g1 - g2)
+
+cg1 = Zygote.gradient(constr, x0)[1]
+cg2 = FiniteDifferences.grad(central_fdm(5, 1), constr, x0)[1]
+norm(cg1 - cg2)
 
 m = Model(obj) # create optimization model
 addvar!(m, zeros(length(x0)), ones(length(x0))) # setup optimization variables
 Nonconvex.add_ineq_constraint!(m, constr) # setup volume inequality constraint
+
 @time r = Nonconvex.optimize(m, TOBSAlg(), x0; options = TOBSOptions())
+# @time r = Nonconvex.optimize(m, IpoptAlg(), x0; options = IpoptOptions())
+# @time r = Nonconvex.optimize(m, MMA87(), x0; options = MMAOptions())
diff --git a/test/runtests.jl b/test/runtests.jl
index 305134f..a7862ea 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,3 +21,7 @@ using ReliabilityOptimization, Test, FiniteDifferences, Zygote
         @test norm(g1 - g2) < 1e-7
     end
 end
+
+@testset "Multi-material TopOpt example" begin
+    include("example.jl")
+end