Open-Sn / opensn / 18928565313

FieldFunctionInterpolationPoint::Create()

  auto ffi = std::make_shared<FieldFunctionInterpolationPoint>();

  field_func_interpolation_stack.emplace_back(ffi);

  return ffi;

}

FieldFunctionInterpolationPoint::Initialize()

  log.Log0Verbose1() << "Initializing point interpolator.";

  if (field_functions_.empty())

    throw std::logic_error("Unassigned field function in point field function interpolator.");

  const auto& grid = field_functions_.front()->GetSpatialDiscretization().GetGrid();

  std::vector<uint64_t> cells_potentially_owning_point;

  for (const auto& cell : grid->local_cells)

    const auto& vcc = cell.centroid;

    const auto& poi = point_of_interest_;

    const auto nudged_point = poi + 1.0e-6 * (vcc - poi);

    if (grid->CheckPointInsideCell(cell, nudged_point))

      cells_potentially_owning_point.push_back(cell.global_id);

  std::vector<uint64_t> recvbuf;

  mpi_comm.all_gather(cells_potentially_owning_point, recvbuf);

  if (recvbuf.empty())

    throw std::logic_error("FieldFunctionInterpolationPoint::Initialize: No cell identified "

                           "containing the specified point.");

  uint64_t owning_cell_gid = recvbuf.front();

  for (const uint64_t gid : recvbuf)

    owning_cell_gid = std::min(owning_cell_gid, gid);

  locally_owned_ = false;

  for (const uint64_t gid : cells_potentially_owning_point)

    if (gid == owning_cell_gid)

      locally_owned_ = true;

      owning_cell_gid_ = owning_cell_gid;

      break;

}

FieldFunctionInterpolationPoint::Execute()

  if (not locally_owned_)

    return;

  const auto& ref_ff = *field_functions_.front();

  const auto& sdm = ref_ff.GetSpatialDiscretization();

  const auto& grid = sdm.GetGrid();

  const auto& uk_man = ref_ff.GetUnknownManager();

  const auto uid = 0;

  const auto cid = ref_component_;

  const auto field_data = ref_ff.GetGhostedFieldVector();

  const auto& cell = grid->cells[owning_cell_gid_];

  const auto& cell_mapping = sdm.GetCellMapping(cell);

  const size_t num_nodes = cell_mapping.GetNumNodes();

  std::vector<double> node_dof_values(num_nodes, 0.0);

  for (size_t i = 0; i < num_nodes; ++i)

    const auto imap = sdm.MapDOFLocal(cell, i, uk_man, uid, cid);

    node_dof_values[i] = field_data[imap];

  Vector<double> shape_values(num_nodes, 0.0);

  cell_mapping.ShapeValues(point_of_interest_, shape_values);

  point_value_ = 0.0;

  for (size_t i = 0; i < num_nodes; ++i)

    point_value_ += node_dof_values[i] * shape_values(i);

}

FieldFunctionInterpolationPoint::GetPointValue() const

  double global_point_value = 0.0;

  mpi_comm.all_reduce(point_value_, global_point_value, mpi::op::sum<double>());

  return global_point_value;