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++)
    {
        phi_i->vector()->setitem(i,1);
        for(int j=0; j < size; j++)
        {
            phi_j->vector()->setitem(i,1);
            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); 
                x[0]=points_coord_x[k];
                x[1]=points_coord_y[k];
                x[2]=points_coord_z[k];
                dolfin::Array<double> vali(1);
                dolfin::Array<double> valj(1);
                phi_i->eval(vali,x);
                phi_j->eval(valj,x);
                val += vali[0]*valj[0];
            }
            D.setitem(ij,val);
            val = 0;
            phi_j->vector()->zero();
        }
        phi_i->vector()->zero();
    }
}

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();
    tree->build(*mesh,3);
    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;
            cell.get_coordinate_dofs(coordinate_dofs);
            std::vector<double> basis(1);
            ufc::cell ufc_cell;
            cell.get_cell_data(ufc_cell);
            cell.get_cell_topology(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, basis.data(), pp.coordinates(),
                                        coordinate_dofs.data(),
                                        ufc_cell.orientation);

                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.