Assigning nodal values to a function only at the boundary nodes

I have been struggling about having a list of node indices:

nodes = [1, 10, 22, ..., 32]

of nodes at the boundary and an associated list of values

node_values = [10.2, 1.2,0.1, ..., -0.2]

and how to assign this to a function defined on the mesh, in order to use it as boundary conditions, instead of using an analytical function or a constant.

After searchin thoroughly I manage to find a solution that works for low order basis functions, for both scalars and vector. I put here my solution to see if this is the way to do it and how can I extend it to higher order basis functions. Also, if it helps someone:

import numpy
from time import time
from dolfin import *

def dofs_to_verts(mesh,boundaryFunction,boundaryValue):
    """
    Compute the mapping from vertices to degrees of freedom only at the boundary.
    This saves time, since we do not have to loop over all the cells of the domain.
    If the for loop was replaced by a simple c++ function call, it would be much faster.

    INPUTS ::
        mesh :: The mesh where we want to compute the mapping between vertices
                and degrees of freedom.

        boundaryFunction :: A MeshFunction which assigns integer values to
                            edges of the mesh.

        boundaryValue :: An integer value which we use to identify the boundary
                         edges we want to use. The MeshFunction should have
                         this value assigned to some of the edges.
    """

    # get the connectivity of the mesh
    mesh.init(1) # initialize edges
    mesh_topology = mesh.topology() # get the mesh topology
    conn12 = mesh_topology(1,2) # get the cells to which an edge belongs

    # get the number of edges
    n_edges = boundaryFunction.size()

    # get the indices of the boundary edges, as given by the boundaryFunction
    boundaryEdges = numpy.arange(n_edges)
    boundaryEdges = boundaryEdges[boundaryFunction.array() == boundaryValue]

    # Define the function spaces associated to the mesh
    # for now only CG of order 1 works for sure
    V = FunctionSpace(mesh, "CG", 1)
    Vv = VectorFunctionSpace(mesh, "CG", 1, dim=2)

    # Allocate memory space of the array that stores the mapping from
    # vertex number to degree of freedom number for both scalar and vector
    # valued functions (note that for vector valued functions we need to indices
    # since the x component and the y component have different indices)
    vert_to_dofs = numpy.zeros(mesh.num_vertices(), dtype=numpy.uintp)
    vectordofs_to_vert = numpy.zeros([2,mesh.num_vertices()], dtype=numpy.uintp)

    # Get the degree of freedom maps
    dm = V.dofmap()
    dms = [Vv.sub(i).dofmap() for i in range(2)] # the i=0 refers to x component
                                                 # the i=1 refers to y component

    # get the cells of the mesh
    mesh_cells = mesh.cells()

    # Since we are only interested in the boundary vertices we just loop over
    # the boundary edges as given by the MeshFunction and stored in boundaryEdges
    # each time we loop we get the numbering of the degrees of freedom of the
    # vertices of the cell that contains the boundary edge and we construct
    # the mapping from vertex number to degree of freedom number
    for edge in boundaryEdges:
        # get the cell index to which the edge belongs
        cell_ind = conn12(edge)[0]
        # get the vertices of the cell
        vert_inds = mesh_cells[cell_ind]

        # add the degree of freedom numbering to the vertices
        vert_to_dofs[vert_inds] = numpy.array(dm.cell_dofs(cell_ind),dtype=numpy.uintp)

        # now we do the same but for vector valued quantites
        # since, in this case, our vector have two quantities (the x and y components)
        # we have to iterate over the two sub dof maps
        for i, (dms_i, dmcs_i) in enumerate(zip(dms, dms)):
            vectordofs_to_vert[i,vert_inds] = numpy.array(dms_i.cell_dofs(cell_ind),dtype=numpy.uintp)

    # Return vertex to dof mapping
    return vert_to_dofs, vectordofs_to_vert

# this is an example of how to assign values at the nodes for both a scalar
# and a vector valued function to a function

# create a unit square mesh (any other mesh works)
mesh = UnitSquareMesh(20,20)

# --------------------------------------------------------------------------
# this part is optional, if you already know the indices of the vertices at
# the boundary just skip this part

# define the boundary mesh
boundaryMesh = BoundaryMesh(mesh)

# if you already don't have the boundary nodes indices
# get the boundary nodes indices
nodes = boundaryMesh.vertex_map().array()

# --------------------------------------------------------------------------

# --------------------------------------------------------------------------
# This part is just to generate a boundary function which assigns a value
# boundary_value to the boundary edges

boundary_value = 10
# get the global number of the boundary edges
edge_map = boundaryMesh.cell_map().array()
# get the number of edges at the boundary
n_boundaryEdges = edge_map.size

# define the boundary function which defines the boundary
boundaryFunction = MeshFunction('size_t',mesh,1)
n_edges = boundaryFunction.array().size # the total number of edges

# allocate memory space to boundary values
boundary_values = numpy.zeros(n_edges,dtype='uintp')
# assign the boundary_value to the array of boundary_values
boundary_values[edge_map] = boundary_value
# assign the boundary_values to the boundary_function
boundaryFunction.set_values(boundary_values)

# --------------------------------------------------------------------------

# generate some random values for the values to assign to the scalar function at the nodes
node_values_scalar = numpy.random.rand(nodes.size)
# generate some random values for the values to assign to the vector function at the nodes
node_values_vector = numpy.random.rand(2,nodes.size)

start = time()
# get the mapping of vertices at boundary to degrees of freedom
vert_to_dofs, vectordofs_to_vert = dofs_to_verts(mesh,boundaryFunction,boundary_value)

print "It took :: %fs" %(time()-start)

# assign vertex values at the boundary for the scalar function

# generate the function space
V = FunctionSpace(mesh,'Lagrange',1)
# generate the scalar function
myFunction = Function(V)
# get the vector containing the data values at the degrees of freedom
myFunction_vector = myFunction.vector()
# use the mapping from vertices (nodes) to degrees of freedom to assign the
# vertex values of the scalar function
myFunction_vector[vert_to_dofs[nodes]] = node_values_scalar

plot(myFunction)

# assign vertex values at the boundary for the vector function

# generate the function space
Vv = VectorFunctionSpace(mesh, "Lagrange", 1, dim=2)
# generate the vector function
myVectorFunction = Function(Vv)
# get the vector containing the data values at the degrees of freedom
myVectorFunction_vector = myVectorFunction.vector()
# use the mapping from vertices (nodes) to degrees of freedom to assign the
# vertex values of the scalar function
myVectorFunction_vector[vectordofs_to_vert[0,nodes]] = node_values_vector[0,:] # x component
myVectorFunction_vector[vectordofs_to_vert[1,nodes]] = node_values_vector[1,:] # y component

plot(myVectorFunction)

interactive()