Incorporating comments from @RonRahaman. Refs enrico-dev#17

include/enrico/openmc_nek_driver.h

@@ -125,11 +125,6 @@ class OpenmcNekDriver : public CoupledDriver {
  //! ordered according to an MPI_Gatherv operation on Nek5000's local elements.
  xt::xtensor<double, 1> elem_volumes_;
 
- //! The dimensionless densities of Nek's global elements
- //! These are **not** ordered by Nek's global element indices. Rather, these are
- //! ordered according to an MPI_Gatherv operation on Nek5000's local elements.
- xt::xtensor<double, 1> elem_densities_;
-
  //! Map that gives a list of Nek element global indices for a given OpenMC material
  //! index. The Nek global element indices refer to indices defined by the MPI_Gatherv
  //! operation, and do not reflect Nek's internal global element indexing.

src/openmc_nek_driver.cpp

@@ -265,15 +265,20 @@ void OpenmcNekDriver::init_densities()
  // each Nek element is fully contained within an OpenMC cell, i.e. Nek
  // elements are not split between multiple OpenMC cells.
  for (int i = 0; i < openmc_driver_->cells_.size(); ++i) {
- if (cell_fluid_mask_[i] == 1) {
- auto& c = openmc_driver_->cells_[i];
- const auto& global_elems = mat_to_elems_.at(c.material_index_);
+ auto& c = openmc_driver_->cells_[i];
+ const auto& global_elems = mat_to_elems_.at(c.material_index_);
 
+ if (cell_fluid_mask_[i] == 1) {
  for (int elem: global_elems) {
  double rho = c.get_density();
  densities_[elem] = rho;
  densities_prev_[elem] = rho;
  }
+ } else {
+ for (int elem: global_elems) {
+ densities_[elem] = 0.0;
+ densities_prev_[elem] = 0.0;
+ }
  }
  }
  }