Output statistics about surface facets rather than vertices of facets

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))

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