 # 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=points_coord_x[k];
x=points_coord_y[k];
x=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*valj;
}
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,localPoint,localPoint);
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;
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;
}

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.