In order to create the data file for a secondary simulation (Monte Carlo), I would like to output information about the surface facets of my 3D mesh into a txt file. Specifically, I would like to create a four column file with the following:

- center of mass coordinate (x,y,z) of the surface element, rather than the coordinates of its corresponding vertices.
- normal vector coordinates of the surface element itself, rather than the normal to the Facet edges that is output from using
`FacetNormal`

- average value of my solution,
`T`

on the given surface element. - area of the corresponding surface element

I have seen that there is an inbuilt function called `FacetArea`

but I am unsure how this can be used to list the corresponding values of the surface area. How can I accomplish this?

The code I currently have to output the spatial coordinates, values of `T`

and the normal vectors of the vertices as output from `FacetNormal(mesh)`

is given below:

```
from dolfin import *
import numpy as np
from mpi4py import MPI
mesh = UnitCubeMesh(1,1,1)
Space = FunctionSpace(mesh,'P', 1)
Temperature = project(Expression("x[0]+0.1*x[1] + 0.1*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)))
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u,v)*ds
l = inner(n, v)*ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)
num_dofs_per_component = int(V.dim()/V.num_sub_spaces())
num_sub_spaces = V.num_sub_spaces()
vector = np.zeros((num_sub_spaces, num_dofs_per_component))
A.ident_zeros()
nh = Function(V)
solve(A, nh.vector(), L)
x = V.sub(0).collapse().tabulate_dof_coordinates()
for i in range(num_sub_spaces):
vector[i] = nh.sub(i, deepcopy=True).vector().get_local()
vector = vector.T
for coord, vec, temp in zip(x_Space, vector, T_values):
print(coord, vec, temp)
```

Thank you in advance