NearestPointBase.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

// MOOSE includes
#include "SpatialUserObjectFunctor.h"
#include "Enumerate.h"
#include "DelimitedFileReader.h"
#include "LayeredBase.h"
#include "FEProblemBase.h"
#include "Positions.h"

// Forward Declarations
class UserObject;

template <typename UserObjectType, typename BaseType>
class NearestPointBase : public SpatialUserObjectFunctor<BaseType>
{
public:
  static InputParameters validParams();

  NearestPointBase(const InputParameters & parameters);
  ~NearestPointBase();

  virtual void initialize() override;
  virtual void execute() override;
  virtual void finalize() override;
  virtual void threadJoin(const UserObject & y) override;

  virtual Real spatialValue(const Point & p) const override;

  virtual const std::vector<Point> & getPoints() const { return _points; }

  virtual const std::vector<Point> spatialPoints() const override;

protected:
  void fillPoints();

  UserObjectType & nearestUserObject(const Point & p) const;

  std::vector<Point> _points;
  std::vector<std::unique_ptr<UserObjectType>> _user_objects;

  // To specify whether the distance is defined based on point or radius
  const unsigned int _dist_norm;
  // The axis around which the radius is determined
  const unsigned int _axis;

  using BaseType::_communicator;
  using BaseType::_current_elem;
  using BaseType::isParamValid;
  using BaseType::processor_id;
  using MooseBase::name;
};

template <typename UserObjectType, typename BaseType>
InputParameters
NearestPointBase<UserObjectType, BaseType>::validParams()
{
  InputParameters params = SpatialUserObjectFunctor<BaseType>::validParams();

  params.addParam<std::vector<Point>>("points",
                                      "Computations will be lumped into values at these points.");
  params.addParam<FileName>("points_file",
                            "A filename that should be looked in for points. Each "
                            "set of 3 values in that file will represent a Point. "
                            "This and 'points' cannot be both supplied.");
  params.addParam<PositionsName>(
      "positions_object",
      "The name of a Positions object that will contain "
      "the locations. This, 'points' and 'points(_file)' cannot be both supplied. "
      "Note that only the vector of initial Positions are used at this time. "
      "Updates to the 'positions' vector are not supported.");

  MooseEnum distnorm("point=0 radius=1", "point");
  params.addParam<MooseEnum>(
      "dist_norm", distnorm, "To specify whether the distance is defined based on point or radius");
  MooseEnum axis("x=0 y=1 z=2", "z");
  params.addParam<MooseEnum>("axis", axis, "The axis around which the radius is determined");

  params.addParamNamesToGroup("points points_file dist_norm axis", "Points and distance to points");

  // Add in the valid parameters
  params += UserObjectType::validParams();

  return params;
}

template <typename UserObjectType, typename BaseType>
NearestPointBase<UserObjectType, BaseType>::NearestPointBase(const InputParameters & parameters)
  : SpatialUserObjectFunctor<BaseType>(parameters),
    _dist_norm(this->template getParam<MooseEnum>("dist_norm")),
    _axis(this->template getParam<MooseEnum>("axis"))
{
  if (this->template getParam<MooseEnum>("dist_norm") != "radius" &&
      parameters.isParamSetByUser("axis"))
    this->paramError("axis", "'axis' should only be set if 'dist_norm' is set to 'radius'");

  fillPoints();

  _user_objects.resize(_points.size());

  // Build each of the UserObject objects that we will manage manually
  // If you're looking at this in the future and want to replace this behavior,
  // _please_ don't do it. MOOSE should manage these objects.
  for (MooseIndex(_points) i = 0; i < _points.size(); ++i)
  {
    const auto uo_type = MooseUtils::prettyCppType<UserObjectType>();
    auto sub_params = this->_app.getFactory().getValidParams(uo_type);
    sub_params.applyParameters(parameters, {}, true);

    const auto sub_name = name() + "_sub" + std::to_string(i);
    auto uo = this->_app.getFactory().template createUnique<UserObjectType>(
        uo_type, sub_name, sub_params, this->_tid);
    _user_objects[i] = std::move(uo);
  }
}

template <typename UserObjectType, typename BaseType>
NearestPointBase<UserObjectType, BaseType>::~NearestPointBase()
{
}

template <typename UserObjectType, typename BaseType>
void
NearestPointBase<UserObjectType, BaseType>::fillPoints()
{
  if (isParamValid("points"))
  {
    _points = this->template getParam<std::vector<Point>>("points");
    if (isParamValid("points_file"))
      this->paramError("points_file", "Cannot be specified together with 'points'");
    if (isParamValid("positions_object"))
      this->paramError("positions_object", "Cannot be specified together with 'points'");
  }
  else if (isParamValid("points_file"))
  {
    const FileName & points_file = this->template getParam<FileName>("points_file");

    MooseUtils::DelimitedFileReader file(points_file, &_communicator);
    file.setFormatFlag(MooseUtils::DelimitedFileReader::FormatFlag::ROWS);
    file.read();
    _points = file.getDataAsPoints();

    if (isParamValid("positions_object"))
      this->paramError("positions_object", "Cannot be specified together with 'points_file'");
  }
  else if (isParamValid("positions_object"))
  {
    const auto & positions_name = this->template getParam<PositionsName>("positions_object");
    const auto problem = this->template getCheckedPointerParam<FEProblemBase *>("_fe_problem_base");
    _points = problem->getPositionsObject(positions_name).getPositions(/*initial_positions*/ true);
  }
  else
    mooseError(
        name(),
        ": You need to supply either 'points', 'points_file' or 'positions_object' parameter.");
}

template <typename UserObjectType, typename BaseType>
void
NearestPointBase<UserObjectType, BaseType>::initialize()
{
  for (auto & user_object : _user_objects)
    user_object->initialize();
}

template <typename UserObjectType, typename BaseType>
void
NearestPointBase<UserObjectType, BaseType>::execute()
{
  nearestUserObject(_current_elem->vertex_average()).execute();
}

template <typename UserObjectType, typename BaseType>
void
NearestPointBase<UserObjectType, BaseType>::finalize()
{
  for (auto & user_object : _user_objects)
    user_object->finalize();
}

template <typename UserObjectType, typename BaseType>
void
NearestPointBase<UserObjectType, BaseType>::threadJoin(const UserObject & y)
{
  auto & npla = static_cast<const NearestPointBase &>(y);

  for (MooseIndex(_user_objects) i = 0; i < _user_objects.size(); ++i)
    _user_objects[i]->threadJoin(*npla._user_objects[i]);
}

template <typename UserObjectType, typename BaseType>
Real
NearestPointBase<UserObjectType, BaseType>::spatialValue(const Point & p) const
{
  return nearestUserObject(p).spatialValue(p);
}

template <typename UserObjectType, typename BaseType>
UserObjectType &
NearestPointBase<UserObjectType, BaseType>::nearestUserObject(const Point & p) const
{
  unsigned int closest = 0;
  Real closest_distance = std::numeric_limits<Real>::max();

  for (auto it : Moose::enumerate(_points))
  {
    const auto & current_point = it.value();

    Real current_distance;
    if (_dist_norm == 0)
      // the distance is computed using standard norm
      current_distance = (p - current_point).norm();
    else
    {
      // the distance is to be computed based on radii
      // in that case, we need to determine the 2 coordinate indices
      // that define the radius
      unsigned int i = 0;
      unsigned int j = 1;

      if (_axis == 0)
        i = 2;
      else if (_axis == 1)
        j = 2;

      current_distance = std::abs(
          std::sqrt(p(i) * p(i) + p(j) * p(j)) -
          std::sqrt(current_point(i) * current_point(i) + current_point(j) * current_point(j)));
    }

    if (current_distance < closest_distance)
    {
      closest_distance = current_distance;
      closest = it.index();
    }
  }

  return *_user_objects[closest];
}

template <typename UserObjectType, typename BaseType>
const std::vector<Point>
NearestPointBase<UserObjectType, BaseType>::spatialPoints() const
{
  std::vector<Point> points;

  for (MooseIndex(_points) i = 0; i < _points.size(); ++i)
  {
    auto layered_base = dynamic_cast<LayeredBase *>(_user_objects[i].get());
    if (layered_base)
    {
      const auto & layers = layered_base->getLayerCenters();
      auto direction = layered_base->direction();

      for (const auto & l : layers)
      {
        Point pt = _points[i];
        pt(direction) = l;
        points.push_back(pt);
      }
    }
  }

  return points;
}