I would start with considering:
which you can restrict the facets in question with:
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(2, 3)
V = FunctionSpace(mesh, "CG", 1)
u = project(Expression("x[0]+3*x[1]", degree=1), V)
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
mf = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
left = Left()
left.mark(mf, 1)
facets = mf.where_equal(1)
send_coordinates = []
send_values = []
send_global_idx = []
global_vertices = mesh.topology().global_indices(0)
vertex_values = u.compute_vertex_values()
coordinates = mesh.coordinates()
for facet_index in facets:
facet = Facet(mesh, facet_index)
vertices = facet.entities(0)
for vertex in vertices:
send_global_idx.append(global_vertices[vertex])
send_coordinates.append(coordinates[vertex])
send_values.append(vertex_values[vertex])
send_coordinates = np.asarray(
send_coordinates).reshape(-1, coordinates.shape[1])
gathered_coordinates = mesh.mpi_comm().gather(send_coordinates, root=0)
gathered_values = mesh.mpi_comm().gather(send_values, root=0)
global_indices = mesh.mpi_comm().gather(send_global_idx, root=0)
if mesh.mpi_comm().rank == 0:
unique_vertices = np.sort(np.unique(np.hstack(global_indices)))
facet_values = np.zeros(len(unique_vertices))
facet_coordinates = np.zeros((len(unique_vertices), coordinates.shape[1]))
c = np.vstack(gathered_coordinates)
v = np.hstack(gathered_values)
i = np.hstack(global_indices)
for coord, value, index in zip(c, v, i):
local_idx = np.flatnonzero(unique_vertices == index)[0]
facet_values[local_idx] = value
facet_coordinates[local_idx] = coord
for (coord, value) in zip(facet_coordinates, facet_values):
print(coord, value)
yielding
mpirun -n 4 python3 script.py
[1. 0.] 0.9999999999999999
[1. 0.33333333] 2.0000000000000004
[1. 0.66666667] 3.000000000000001
[1. 1.] 4.0000000000000036
You can easily store this to a text file using the normal file handler in Python