How to loop over facets on boundaries and their coodinates to export .pos file to Gmsh ?

Question

Dear All,

I am currently trying to solve a 3D Laplace equation for the scalar potenial u over 2 concentric cubes (w/ 2.6 and 20 units in length) meshed in Gmsh and converted using dolfin-convert:
- Outer cube surface (\Gamma_{1}) is set at potential u1=0 (Dirichlet BC)
- inner cube bottom (\Gamma_{2}) and top (\Gamma_{3}) surfaces are respectively set at constant Neumann BC 1 and -1. All other inner surfaces (\Gamma_{4}) to (\Gamma_{7}) have 0 Neumann BC.
I (try to) solve the problem and would like to export the solution to Gmsh in .pos format where
for each triangle with corners (V1,V2,V3) in the mesh boundary one has a line like:
ST(V1x,V1y,V1z,V2x,V2y,V2z,V3x,V3y,V3z){u1,u2,u3};
with ui the potential at corner Vi of coordinates (Vix,Viy,Viz).
To write such a file, I need to loop over all triangles on the 6 inner boundaries and for each triangle I need to retrieve the corners coordinates and associated potential solution, but I can't figure out how to access these values.
The Python code used is:

from dolfin import *

Create mesh and define function space

mesh = Mesh("single_layer_coil.xml")

subdomains = MeshFunction("size_t", mesh, "single_layer_coil_physical_region.xml")

boundaries = MeshFunction("size_t", mesh, "single_layer_coil_facet_region.xml")

V = FunctionSpace(mesh, "CG", 1)

Define Neumann boundary conditions

Mark bottom (-Z) boundary facets as subdomain 2

class BottomNeumannBoundary(SubDomain):

    def inside(self, x, on_boundary):

        tol = 1E-14   # tolerance for coordinate comparisons

        return on_boundary and abs(x[2] + 1.3) < tol and x[0] <= 1.3 and x[0] >= -1.3 and x[1] <= 1.3 and x[1] >= 1.3

Gamma_N2 = BottomNeumannBoundary()

Gamma_N2.mark(boundaries, 2)

g2 = Constant(1)

Mark top (+Z) boundary facets as subdomain 3

class TopNeumannBoundary(SubDomain):

    def inside(self, x, on_boundary):

        tol = 1E-14   # tolerance for coordinate comparisons

        return on_boundary and abs(x[2] - 1.3) < tol and x[0] <= 1.3 and x[0] >= -1.3 and x[1] <= 1.3 and x[1] >= -1.3

Gamma_N3 = TopNeumannBoundary()

Gamma_N3.mark(boundaries, 3)

g3 = Constant(-1)

Define Dirichlet boundary condition values

Physical Surface(0)

u1  = Constant(0)

bc1 = DirichletBC(V, u1, boundaries, 1)

Define surface element

ds = Measure('ds')[boundaries]

Define variational problem

u = TrialFunction(V)

v = TestFunction(V)

f = Constant(0.0)

a = dot(grad(u), grad(v))dx

L = fvdx - g2vds(2) - g3v*ds(3)

Compute solution

u = Function(V)

solve(a == L, u, [bc1])

Write solution to file

File("u.pvd") << u

Plot solution

plot(u, interactive=True)

Extract solution over mesh vertices

u_nodal_values = u.vector()

u_array = u_nodal_values.array()

coor = mesh.coordinates()

if mesh.num_vertices() == len(u_array):

for i in range(mesh.num_vertices()):

print 'u(% 8g,% 8g,% 8g) = %g' % (coor[i][0], coor[i][1], coor[i][2], u_array[i])

Another related question: how can one plot the solution over some parts of the mesh ?
Thanks a lot for any help !