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.
I faced the same problem and thanks for sharing your solution. However, something confused me. What is the difference between coordinates and coordinates_dofs? Is the coordinates reshaped coordinates_dof?
Thank you very much!
The final code is a bit different than the one posted above:
void point_cloud_domain::populate_delta_mass_matrix(std::shared_ptr<dolfin::PETScMatrix> A,
std::shared_ptr<dolfin::Mesh> mesh, double beta, double initalVol) const
{
auto Vh = std::make_shared<varfs::FunctionSpace>(mesh);
auto phi = std::make_shared<dolfin::Function>(Vh);
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);
dolfin::Point pp(localPoint[0],localPoint[1],localPoint[2]);
dolfin::Array<double> x(3); x[0] = pp[0]; x[1] = pp[1]; x[2] = pp[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);
double ratioVol = initalVol/cell.volume();
std::vector<long unsigned int> indeces = ufc_cell.entity_indices[0];
dolfin::Array<double> basis_data(1);
std::vector<double> valueAtPoint(4);
for (std::size_t i = 0; i < 4; ++i)
{
phi->vector()->zero();
phi->vector()->setitem(indeces[i],1.0);
phi->eval(basis_data,x);
valueAtPoint[i] = basis_data[0];
}
for(int i=0; i < 4; i++)
{
for(int j=0; j < 4; j++)
{
double value = ratioVol*beta*valueAtPoint[i]*valueAtPoint[j];
if(value > 0)
{
std::pair<dolfin::la_index,dolfin::la_index> ij(indeces[i],indeces[j]);
A->add_local(&value,1,&ij.first,1,&ij.second);
}
}
}
}
}
A->apply("add");
std::cout << "Delta mass matrix populated" << std::endl;
}
This code works. This is because I honestly don’t know what evaluate_basis function does in practice, and it produced nonsense results.
I don’t remember why I needed coordinate_dofs, but it seems not to be useful anymore. The coordinates() method returns the list of the coordinates of the vertices of a mesh such that coord[i*3+j] is the jth component of the ith vertex (for 3D meshes).
Maybe this code can be optimised far more, but without any C++ documentation I think this is the best result I can get.
1 Like