But I need the coordinates of a part of the boundary only, subdomain 1, say, and it is difficult in this case to sort these out afterwards. So, how can I get the coordinates of just the subboundary directly?

By the way, I also would like to evaluate a function u exactly at these nodes. So, how can I relate the entries of u.array()[:] with those boundary nodes?

That is not an error, that is simply from the fact that a CG-2 space has degrees of freedom on facets as well as vertices. Therefore you cannot use vertex_to_dof_map for this operation. You need a more complex operation for this:

from dolfin import *
import numpy as np
class boundary(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 0.5 and on_boundary
mesh = UnitSquareMesh(5,5)
V = FunctionSpace(mesh, "CG", 2)
v = project(Expression("x[0]+x[1]", degree=2),V)
mf = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundary().mark(mf, 1)
mesh.init(2,1)
dofs = []
cell_to_facets = mesh.topology()(2,1)
for cell in cells(mesh):
facets = cell_to_facets(cell.index())
for facet in facets:
if mf[facet] == 1:
dofs_ = V.dofmap().entity_closure_dofs(mesh, 1, [facet])
for dof in dofs_:
dofs.append(dof)

The MeshTopology() class is not adequately documented, and Its not clear what the commands mesh.init(2,1) and mesh.topology()(2,1) do. How we can get more information to use in different applicaions?

This function creates the connectivity between entities of topological dimension 2 (triangles) to entities of topological dimension 1 (lines) using indices local to a process. Each triangle will have three line segments connected, i.e.

gives you a map from each cell (local to process) to the three line segments (index local to process).