Hi,
I’m wondering if there is an easy way to get the map from DOF to elements containing that DOF. It’s easy when using P1 elements since the DOF are the mesh vertices and fenics makes a map from vertices to cells available (via the built-in mesh entity iterators), but I would prefer a more general approach. I know I can construct the map myself by looping through the elements and getting the associated DOF (i.e., using the cell -> DOF map to get the inverse), but I’m wondering if there is another way already built into fenics.
I guess my question is similar to this one, and the solutions I have come up with for both the P1 case and the general case are mentioned in the comments on the first reply. This suggests that maybe there is no better way, but I will leave this open just in case someone has a suggestion.
Update - this was my solution (it’s C++ code):
auto mesh = std::make_shared<UnitSquareMesh>(100, 100);
auto V = std::make_shared<poisson::Form_J_FunctionSpace_2>(mesh);
auto dofmap = (*V).dofmap();
std::vector<std::vector<size_t>> cell_data((*V).dim());
for(CellIterator eit(*mesh); !eit.end(); ++eit){ //loop over cells
size_t index = (*eit).index();
Eigen::Map<const Eigen::Array<int, -1, 1>> dof = (*dofmap).cell_dofs(index); //get indices of dofs on each cell
size_t val;
for(int i = 0; i<dof.size();i++){// loop over the dofs
val = dof(i);
if(cell_data[val].size() == 0)
cell_data[val].push_back(index);
else {// make sure that the index is not already in the list before pushing it on
for(int j=0; j<cell_data[val].size(); j++){
if(cell_data[val][j] == index)
break;
else if(j + 1 == cell_data[val].size())
cell_data[val].push_back(index);
};
};
};
};
It’s a bit hairy, but I guess it works.