Access getMesh() protected method for non meshgenerators

[lipohlmann]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Hello all!

I am working to develop a SurfaceFluxGenerator in Magpie. To do this, I am having the user specify the boundary and the mesh on which the boundary exists. Within the class, methods have been developed to determine the boundary surface area and a weighted CDF for particle sampling. The issue is that the getMesh() method is not available because this userobject does not inherit from MeshGenerator. Would anyone have any suggestions for how to go about getting access to this and its dependent methods in order to pull in the mesh correctly? Or would there be a better way to access the boundary mesh?

Header file:

/**********************************************************************/
/*                     DO NOT MODIFY THIS HEADER                      */
/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
/*                                                                    */
/*            Copyright 2017 Battelle Energy Alliance, LLC            */
/*                        ALL RIGHTS RESERVED                         */
/**********************************************************************/

#pragma once

#include "PKAGeneratorBase.h"

/**
 * Starts PKAs at a fixed point in a fixed direction
 */
class PKASurfaceFluxGenerator : public PKAGeneratorBase
{
public:
  static InputParameters validParams();

  PKASurfaceFluxGenerator(const InputParameters & parameters);

  virtual void appendPKAs(std::vector<MyTRIM_NS::IonBase> &,
                          const MyTRIMRasterizer::PKAParameters &,
                          const MyTRIMRasterizer::AveragedData &) const;

protected:
  /// provides a mean to override the angular distribution of the PKAs in derived class
  virtual void setDirection(MyTRIM_NS::IonBase & ion) const;

  /// the direction along which the PKAs start
  const RealVectorValue _direction;

  /// Uses point locator to determine the element id of the elemnt _point is in
  virtual void updateCachedElementID();

  /// number of PKAs to be generated per cm2-s
  const unsigned int _flux;

  /// time step
  const Real _dt;

  /// Mesh that comes from another generator
  const MeshBase * _mesh;

  /// boundary name
  const BoundaryName _boundary;

  /// surface area of boundary
  const Real _boundary_surface_area;

  /// the location from which to start PKAs
  Point _point;

  /// PKA nuclear charge
  const unsigned int _Z;

  /// PKA mass
  const Real _m;

  /// PKA Energy (in eV)
  const Real _E;

  /// cumulative probability list with element IDs
  const std::vector<std::pair<Real, dof_id_type>> _prob_elem_pairs;

  /// point locator to determine element pointers form locations
  std::unique_ptr<PointLocatorBase> _pl;

  /// the element id of the element containing _point
  dof_id_type _elem_id;

  Real boundarySurfaceArea(const BoundaryName & boundary, const std::unique_ptr<MeshBase> & mesh);

  std::vector<std::pair<Real, dof_id_type>>
  volumeWeightedElemDist(const BoundaryName & boundary, const std::unique_ptr<MeshBase> & mesh);
};

Source File:

/**********************************************************************/
/*                     DO NOT MODIFY THIS HEADER                      */
/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
/*                                                                    */
/*            Copyright 2017 Battelle Energy Alliance, LLC            */
/*                        ALL RIGHTS RESERVED                         */
/**********************************************************************/

#include "PKASurfaceFluxGenerator.h"
#include "MagpieUtils.h"

registerMooseObject("MagpieApp", PKASurfaceFluxGenerator);

InputParameters
PKASurfaceFluxGenerator::validParams()
{
  InputParameters params = PKAGeneratorBase::validParams();
  params.addClassDescription(
      "This PKAGenerator starts particles from a random point in a boundary in a fixed direction.");
  params.addRequiredParam<RealVectorValue>("direction", "The fixed direction the PKAs move along");
  params.addRequiredParam<MeshGeneratorName>(
      "mesh", "The mesh where the boundary to apply the flux to lives.");
  params.addRequiredParam<BoundaryName>("boundary", "The boundary to apply the flux to.");
  params.addParam<unsigned int>(
      "flux", 1e8, "The number of ions (starting PKAs) to strike the surface in ions/cm2-sec");
  params.addRequiredParam<Real>("dt", "Time step size used in time stepper. Value in cm2!");
  params.addRequiredParam<Real>("boundary_surface_area",
                                "Surface area of boundary. Used to determine number of initial "
                                "PKAs generated per timestep. Value in cm2!");
  params.addRequiredParam<Real>("Z", "PKA nuclear charge");
  params.addRequiredParam<Real>("m", "PKA mass in amu");
  params.addRequiredParam<Real>("E", "PKA energy in eV");
  return params;
}

PKASurfaceFluxGenerator::PKASurfaceFluxGenerator(const InputParameters & parameters)
  : PKAGeneratorBase(parameters),
    _direction(getParam<RealVectorValue>("direction")),
    _flux(getParam<Real>("flux")),
    _dt(getParam<Real>("dt")),
    _mesh(getMesh("mesh")),
    _boundary(getParam<BoundaryName>("boundary")),
    _boundary_surface_area(PKASurfaceFluxGenerator::boundarySurfaceArea(_boundary, _mesh)),
    _Z(getParam<Real>("Z")),
    _m(getParam<Real>("m")),
    _E(getParam<Real>("E")),
    _prob_elem_pairs(PKASurfaceFluxGenerator::volumeWeightedElemDist(_boundary, _mesh))
{
}

void
PKASurfaceFluxGenerator::appendPKAs(std::vector<MyTRIM_NS::IonBase> & ion_list,
                                    const MyTRIMRasterizer::PKAParameters & pka_parameters,
                                    const MyTRIMRasterizer::AveragedData &) const
{
  if (_current_elem->id() != _elem_id)
    return;

  unsigned int num_pka = _flux * _dt * _boundary_surface_area;
  if (pka_parameters._recoil_rate_scaling != 1)
    num_pka = std::floor(pka_parameters._recoil_rate_scaling * num_pka + getRandomReal());

  for (unsigned i = 0; i < num_pka; ++i)
  {
    // each fission event generates a pair of recoils
    MyTRIM_NS::IonBase pka;

    // sample fission fragment masses
    pka._Z = _Z;
    pka._m = _m;
    pka._E = _E;

    // the tag is the element this PKA get registered as upon stopping
    // -1 means the PKA will be ignored
    pka._tag = ionTag(pka_parameters, pka._Z, pka._m);
    ;

    // set stopping criteria
    pka.setEf();

    // sample elements
    Real rnd_num = getRandomReal();
    auto it = std::lower_bound(_prob_elem_pairs.begin(),
                               _prob_elem_pairs.end(),
                               rnd_num,
                               [](const std::pair<Real, dof_id_type> & entry, const Real & value)
                               { return entry.first < value; });

    // Get the selected element ID
    dof_id_type selected_elem_id = it->second;

    // Retrieve the element pointer from mesh
    Elem * rnd_elem = _mesh->elem_ptr(selected_elem_id);

    Point _point = MagpieUtils::randomElementPoint(
        *rnd_elem, getRandomPoint()); // sample random point in element
    pka._pos = _point;

    // set random direction for ion 1 and opposite direction for ion 2
    setDirection(pka);

    // add PKA to list
    ion_list.push_back(pka);
  }
}

void
PKASurfaceFluxGenerator::setDirection(MyTRIM_NS::IonBase & ion) const
{
  ion._dir = _direction;
}

Real
PKASurfaceFluxGenerator::boundarySurfaceArea(const BoundaryName & boundary,
                                             const std::unique_ptr<MeshBase> & mesh)
{
  if (!mesh->is_serial())
    mooseError("findCenterPoint not yet implemented for distributed meshes!");

  BoundaryInfo & mesh_boundary_info = mesh->get_boundary_info();
  boundary_id_type boundary_id = mesh_boundary_info.get_id_by_name(std::string_view(boundary));

  // initialize sums
  double volume_sum = 0;

  // loop over all elements in mesh
  for (const auto & elem : mesh->element_ptr_range())
  {
    // loop over all sides in element
    for (const auto side_num : make_range(elem->n_sides()))
    {
      // check if on boundary
      bool on_boundary = mesh_boundary_info.has_boundary_id(elem, side_num, boundary_id);
      if (on_boundary)
        volume_sum += elem->side_ptr(side_num)->volume(); // update running sum
    }
  }
  return volume_sum;
}

std::vector<std::pair<Real, dof_id_type>>
PKASurfaceFluxGenerator::volumeWeightedElemDist(const BoundaryName & boundary,
                                                const std::unique_ptr<MeshBase> & mesh)
{
  if (!mesh->is_serial())
    mooseError("volumeWeightedElemDist not yet implemented for distributed meshes!");

  BoundaryInfo & mesh_boundary_info = mesh->get_boundary_info();
  boundary_id_type boundary_id = mesh_boundary_info.get_id_by_name(std::string_view(boundary));

  double volume_sum = 0.0;
  std::vector<std::pair<dof_id_type, Real>> elem_volumes; // (elem_id, side_volume)
  std::vector<std::pair<Real, dof_id_type>> prob_elem_pairs;

  // Loop over all elements in the mesh
  for (const auto & elem : mesh->element_ptr_range())
  {
    dof_id_type elem_id = elem->id();

    for (const auto side_num : make_range(elem->n_sides()))
    {
      if (mesh_boundary_info.has_boundary_id(elem, side_num, boundary_id))
      {
        Real side_vol = elem->side_ptr(side_num)->volume();
        volume_sum += side_vol;
        elem_volumes.emplace_back(elem_id, side_vol);
      }
    }
  }

  // Normalize into cumulative distribution
  Real cumulative_sum = 0.0;
  for (const auto & [elem_id, side_vol] : elem_volumes)
  {
    cumulative_sum += side_vol / volume_sum;
    prob_elem_pairs.emplace_back(cumulative_sum, elem_id);
  }

  mooseAssert(std::abs(cumulative_sum - 1.0) < libMesh::TOLERANCE,
              "Volumes are not normalized to sum to 1.0!");
  mooseAssert(std::abs(prob_elem_pairs.back().first - 1.0) < libMesh::TOLERANCE,
              "Element probabilities are not normalized to sum to 1.0!");

  return prob_elem_pairs;
}

[GiudGiud]

Hi Liam

Are you using multiple meshes?

If not, then use the _problem or _fe_problem attribute of your object and get the mesh using problem.mesh()

If so, then you can use the "save_in" option of mesh generators to save a mesh, and then retrieve it using the MeshGeneratorSystem::getSavedMesh
The MGSystem can be retrieved from the MooseApp _app attribute of most objects with _app.getMeshGeneratorSystem()

[lipohlmann]

Sorry for the late reply. _fe_problem.mesh() worked well! I'll post the constructor and relevant snippet of the header for anyone who sees this thread:

Constructor in source file:

PKASurfaceFluxGenerator::PKASurfaceFluxGenerator(const InputParameters & parameters)
  : PKAGeneratorBase(parameters),
    _direction(getParam<RealVectorValue>("direction")),
    _flux(getParam<unsigned int>("flux")),
    _dt(getParam<Real>("dt")),
    _mesh(_fe_problem.mesh()),
    _boundary(getParam<BoundaryName>("boundary")),
    _boundary_surface_area(PKASurfaceFluxGenerator::boundarySurfaceArea(_boundary)),
    _Z(getParam<Real>("Z")),
    _m(getParam<Real>("m")),
    _E(getParam<Real>("E")),
    _prob_elem_pairs(PKASurfaceFluxGenerator::volumeWeightedElemDist(_boundary)),
    _pl(_mesh.getPointLocator())
{
  _pl->enable_out_of_mesh_mode();
  updateCachedElementID();
}

From protected attributes of the header:

...
/// Mesh that comes from another generator
MooseMesh & _mesh;

/// boundary name
const BoundaryName _boundary;
...

The fact that mesh() returns a reference to a MooseMesh instead of a pointer through me for a loop at first...