diff --git a/docs/src/examples/methods.jl b/docs/src/examples/methods.jl
index a53a8da..b9aa69f 100644
--- a/docs/src/examples/methods.jl
+++ b/docs/src/examples/methods.jl
@@ -11,21 +11,22 @@ import FileIO # hide
 domain = (0,1,0,1)
 cells = (10,10)
 mesh = GT.cartesian_mesh(domain,cells)
+D = GT.num_dims(mesh)
 GT.label_boundary_faces!(mesh;physical_name="Γd")
 Ω = GT.interior(mesh)
 Γd = GT.boundary(mesh;physical_names=["Γd"])
-n_Γd = GT.unit_normal(Γd,Ω)
+n_Γd = GT.unit_normal(mesh,D-1)
 u = GT.analytical_field(x->sum(x),Ω)
 #TODO these two definitions should be enough
-#∇ = ForwardDiff.gradient
-#Δ = (u,x) -> tr(ForwardDiff.jacobian(y->∇(u,y),x))
-gradient(u) = x->ForwardDiff.gradient(u,x)
-jacobian(u) = x->ForwardDiff.jacobian(u,x)
-laplacian(u) = x-> tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x))
-Δ(u) = GT.call(laplacian,u)
-∇(u) = GT.call(gradient,u)
-Δ(u,x) = Δ(u)(x)
-∇(u,x) = ∇(u)(x)
+∇ = ForwardDiff.gradient
+Δ = (u,x) -> tr(ForwardDiff.jacobian(y->∇(u,y),x))
+#gradient(u) = x->ForwardDiff.gradient(u,x)
+#jacobian(u) = x->ForwardDiff.jacobian(u,x)
+#laplacian(u) = x-> tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x))
+#Δ(u) = GT.call(laplacian,u)
+#∇(u) = GT.call(gradient,u)
+#Δ(u,x) = Δ(u)(x)
+#∇(u,x) = ∇(u)(x)
 tol = 1e-10
 order = 2
 γ = order*(order+1)/10
@@ -54,7 +55,7 @@ el2 = sqrt(sum(GT.∫( x->abs2(eh(x)), dΩ)))
 # ## Nitche Method
 
 V = GT.lagrange_space(Ω,order)
-n_Γd = GT.unit_normal(Γd,Ω)
+n_Γd = GT.unit_normal(mesh,D-1)
 h_Γd = GT.face_diameter_field(Γd)
 dΓd = GT.measure(Γd,2*order)
 a_Γd = (u,v) -> GT.∫(
@@ -90,18 +91,18 @@ el2 = sqrt(sum(GT.∫( x->abs2(eh(x)), dΩ)))
 
 # ## Interior Penalty
 
-mean(u,x) = 0.5*(u(x)[+]+u(x)[-])
-jump(u,n,x) = u(x)[+]*n[+](x) + u(x)[-]*n[-](x)
+mean(u) = 0.5*(u[1]+u[2])
+jump(u,n) = u[2]*n[2] + u[1]*n[1]
 GT.label_interior_faces!(mesh;physical_name="Λ")
 Λ = GT.skeleton(mesh;physical_names=["Λ"])
 dΛ = GT.measure(Λ,2*order)
-n_Λ = GT.unit_normal(Λ,Ω)
+n_Λ = GT.unit_normal(mesh,D-1)
 h_Λ = GT.face_diameter_field(Λ)
 V = GT.lagrange_space(Ω,order;conformity=:L2)
 a_Λ = (u,v) -> GT.∫(
-    x-> (γ/h_Λ(x))*jump(v,n_Λ,x)⋅jump(u,n_Λ,x)-
-    jump(v,n_Λ,x)⋅mean(∇(u),x)-
-    mean(∇(v),x)⋅jump(u,n_Λ,x), dΛ)
+                    x-> (γ/h_Λ(x))*jump(v(x),n_Λ(x))⋅jump(u(x),n_Λ(x))-
+                    jump(v(x),n_Λ(x))⋅mean(∇(u,x))-
+                    mean(∇(v,x))⋅jump(u(x),n_Λ(x)), dΛ)
 a_ip = (u,v) -> a_Ω(u,v) + a_Γd(u,v) + a_Λ(u,v)
 l_ip = l_ni
 p = GT.linear_problem(Float64,V,a_ip,l_ip)
diff --git a/src/symbolics.jl b/src/symbolics.jl
index 67fda28..34967b2 100644
--- a/src/symbolics.jl
+++ b/src/symbolics.jl
@@ -976,6 +976,28 @@ function binary_call_term(d::typeof(ForwardDiff.gradient),a::ComposedWithInverse
     call(\,Jt,gradf)
 end
 
+function binary_call_term(d::typeof(ForwardDiff.jacobian),a::ComposedWithInverseTerm,b::FunctionCallTerm)
+    phi1 = a.arg2.arg1
+    phi2 = b.arg1
+    if phi1 != phi2
+        return binary_call_term_physical_maps(d,a,b,phi1,phi2)
+    end
+    phi = phi1
+    f = a.arg1
+    x = b.arg2
+    J = call(ForwardDiff.jacobian,phi,x)
+    Jf = call(d,f,x)
+    call(/,Jf,J)
+end
+
+#  w(q) = u(ϕ(q))
+#  ∇(w,q) = Jt(ϕ,q)*∇(u,ϕ(q))
+#  Jt(ϕ,q)\∇(w,q) = ∇(u,ϕ(q))
+#
+#
+#  J(w,q) = J(u,ϕ(q))*J(ϕ,q)
+#  J(w,q)/J(ϕ,q) = J(u,ϕ(q))
+
 function binary_call_term_physical_maps(op,a,b,phi1,phi2)
     return BinaryCallTerm(op,a,b)
 end