How to compute the value of a basis function \phi in a specific point?

Good morning!

I go right to the question: I need to compute a mass matrix M such that

M_{ij} = \sum_{k=1}^N \phi_i(x_k) \phi_j(x_k)

Where {x_k} is a set of points (10^4 entries circa). My idea was to directly populate the matrix M, but to do so I need to access each single basis function associated to every mesh element.

I use the C++ API. Unfortunately, the only idea I have is to compute for each vertex a dolfin::Function which has each nodal value set to zero except the one located in i-j. I would like to know wether there are computationally cheaper alternatives because the following code:

void populate_delta_mass_matrix(dolfin::Matrix D,
                     std::shared_ptr<dolfin::Mesh> mesh)
    auto Vh = std::make_shared<varfs::FunctionSpace>(mesh);
    auto phi_i = std::make_shared<dolfin::Function>(Vh); phi_i->vector()->zero();
    auto phi_j = std::make_shared<dolfin::Function>(Vh); phi_j->vector()->zero();
    int size = phi_i->vector()->size();
    for(int i=0; i < size; i++)
        for(int j=0; j < size; j++)
            std::pair<dolfin::la_index,dolfin::la_index> ij(i,j);
            double val = 0;
            for(int k=0; k<size_of_the_point_cloud; k++)
                dolfin::Array<double> x(3); 
                dolfin::Array<double> vali(1);
                dolfin::Array<double> valj(1);
                val += vali[0]*valj[0];
            val = 0;

is too much computationally expensive!

Ok, I found by myself a suitably efficient solution, thanks to the developer comments in Function.cpp:

void populate_delta_mass_matrix(dolfin::Matrix D,
                         std::shared_ptr<dolfin::Mesh> mesh) const
    auto Vh = std::make_shared<varfs::FunctionSpace>(mesh);
    std::shared_ptr<dolfin::BoundingBoxTree> tree = mesh->bounding_box_tree();
    auto element = Vh->element();

    for(int ii=0; ii < this->m_dim; ii++)
        std::vector<double> localPoint = get_point(ii); // it contains the coordinates of the cloud point
        dolfin::Point pp(localPoint[0],localPoint[1],localPoint[2]);
        auto idx = tree->compute_first_entity_collision(pp);
        if(idx > 0 && idx < mesh->num_cells())
            Cell cell(*mesh, idx);
            std::vector<double> coordinate_dofs;
            std::vector<double> basis(1);
            ufc::cell ufc_cell;

            std::vector<long unsigned int> indeces = ufc_cell.entity_indices[0];
            std::vector<double> basis_data(4);

            for (std::size_t i = 0; i < element->space_dimension(); ++i)
                element->evaluate_basis(i,, pp.coordinates(),

                basis_data[i] = basis[0];

            for(int i=0; i < 4; i++)
                for(int j=0; j < 4; j++)
                    std::pair<dolfin::la_index,dolfin::la_index> ij(indeces[i],indeces[j]);
                    double val = D.getitem(ij);
                    D.setitem(ij,val + basis_data[i]*basis_data[j]);

The above function takes in input a mesh and a matrix already populated (entries set to zero). The point cloud belongs to the class of this function, but it is easy to modify it for further purposes.