Actually, one last thing. I have tried stitching together your provided code and the code from @kamensky such that I can output the spatial coordinates of the facet, value of T at that facet and the normal vector components to a txt file, but it does not seem to be working and I am unsure if the spatial coordinate values, T values and normal vector components actually correspond to each other:

```
from dolfin import *
import numpy as np
from FIAT.reference_element import UFCTetrahedron, UFCTriangle
from FIAT import ufc_simplex, Lagrange
mesh = UnitCubeMesh(1, 1, 1)
#mesh = UnitSquareMesh(10, 10)
Space = FunctionSpace(mesh,'P', 1)
Temperature = project(Expression("x[0]+0.1*x[1] -0.4*x[2]", degree=3), Space)
x_Space = Space.tabulate_dof_coordinates().reshape((-1,mesh.geometry().dim()))
T_values = abs(Temperature.vector().get_local().reshape((-1,1)))
ufc_element = UFCTetrahedron() if mesh.ufl_cell() == tetrahedron else UFCTriangle()
tdim = mesh.topology().dim()
fdim = tdim - 1
top = mesh.topology()
mesh.init(tdim, fdim)
mesh.init(fdim, tdim)
f_to_c = top(fdim, tdim)
bndry_facets = []
for facet in facets(mesh):
f_index = facet.index()
cells = f_to_c(f_index)
if len(cells) == 1:
bndry_facets.append(f_index)
# As we are using simplices, we can compute the Jacobian at an arbitrary point
fiat_el = Lagrange(ufc_simplex(tdim), 1)
dphi_ = fiat_el.tabulate(1, np.array([0,]*tdim))
dphi = np.zeros((tdim, Cell(mesh, 0).num_entities(0)), dtype=np.float64)
for i in range(tdim):
t = tuple()
for j in range(tdim):
if i == j:
t += (1, )
else:
t += (0, )
dphi[i] = dphi_[t]
c_to_f = top(tdim, fdim)
n_values = []
for facet in bndry_facets:
cell = f_to_c(facet)[0]
facets = c_to_f(cell)
local_index = np.flatnonzero(facets==facet)[0]
normal = ufc_element.compute_reference_normal(fdim, local_index)
normal_n = normal/np.linalg.norm(normal)
coord_dofs = mesh.cells()[cell]
coords = mesh.coordinates()[coord_dofs]
J = np.dot(dphi, coords).T
Jinv = np.linalg.inv(J)
n_glob = np.dot(Jinv.T, normal_n)
n_glob = n_glob/np.linalg.norm(n_glob)
n_values.append(n_glob)
# print(facet,coords, n_glob)
filename = "output.txt"
data = [x_Space, T_values, n_values]
with open(filename, 'w') as f:
np.savetxt(filename, data, header = 'x,y,z,T(K),nx,ny,nz')
```