Engineering Units

Common Methods

All custom engineering units must implement the following methods (first two methods are described in more detail below):

• expandUnit(): creates the equivalent object using basic CSG components and adds them to the provided CSGBase instance

• getAttributes(): returns a map of necessary unit attributes to their values

• clone(): clones the instance as a unique pointer

All engineering units have additional information retrievable with the following methods defined in the abstract CSGEngUnit class:

• getName(): returns the name of the unit

• getBehavior(): returns what type of object this unit behaves like ("SURFACE", "CELL", or "UNIVERSE")

• getUnitType(): returns the class name as a string for the specific engineering unit

CSGSurfaceEngUnit

To define a surface-like engineering unit, define a class that derives from CSGSurfaceEngUnit. Functionally, these behave like a CSGSurface in that you can use them to define half-spaces when defining CSGRegions for a CSGCell. All methods that are available on CSGSurface objects are also available on CSGSurfaceEngUnit objects except getCoeffs (not applicable to engineering units).

CSGCellEngUnit

To define a cell-like engineering unit, define a class the derives from CSGCellEngUnit. Functionally, these behave like CSGCell objects in that you can add them to a CSGUniverse. Because CSGCellEngUnits are defined by custom human-readable attributes, the standard CSGCell attributes such as a cell fill or region are undefined and unable to be updated via the standard methods available to CSGCell.

CSGUniverseEngUnit

To define a universe-like engineering unit, define a class the derives from CSGUniverseEngUnit. Functionally, these behave like a CSGUniverse object in that you can use them as fills for CSGCells or elements or outer for a CSGLattice. The expectation is that a CSGUniverseEngUnit replaces a collection of CSGCell objects using human-readable information, and therefore, the method addCell is not applicable to this unit type like it is a CSGUniverse.

Example Implementation

Below shows how an engineering unit to define an N-Sided Regular Polygon (CSGNPolygonUnit) as a CSGSurfaceEngUnit is implemented as an example.

#pragma once

#ifdef MOOSE_UNIT_TEST
#include "gtest/gtest.h"
#endif

#include "CSGSurfaceEngUnit.h"

namespace CSG
{

class CSGBase;

/**
 * CSGNPolygonUnit is a CSGSurfaceEngUnit that represents a regular N-sided polygon
 * (prismatic region) with its axis aligned with the z-axis.
 *
 * The polygon is defined by N infinite planes, one per side, each oriented so that
 * its inward normal points toward the center. The interior of the polygon is the
 * intersection of the N negative half-spaces of those planes.
 *
 * Implements:
 *   - expandUnit(): creates N CSGPlane surfaces in _internal_base, one per side
 *   - getExpandedRegion(): returns the intersection of the N half-spaces formed by the planes that
 *     containing the origin
 *   - evaluateSurfaceEquationAtPoint(): returns the maximum plane evaluation over
 *     all N sides; negative means the point is inside, positive means outside
 *   - clone(): returns a deep copy
 *   - getAttributes(): returns n_sides, apothem
 */
class CSGNPolygonUnit : public CSGSurfaceEngUnit
{
public:
  /**
   * @brief Construct a new CSGNPolygonUnit.
   *
   * @param name unique name of the unit
   * @param n_sides number of sides of the regular polygon (must be >= 3)
   * @param apothem distance from the center to the midpoint of each side
   */
  CSGNPolygonUnit(const std::string & name, int n_sides, Real apothem);

  /**
   * @brief Evaluate the polygon's surface equation at the given point.
   *
   * Returns the maximum signed plane-equation value over all N sides. A negative
   * return value means the point lies inside the polygon; positive means outside;
   * zero means on one of the sides.
   *
   * This can be evaluated before expansion using the stored geometric parameters.
   *
   * @param p point to evaluate
   * @return maximum signed distance from the point to any side (negative = inside)
   */
  Real evaluateSurfaceEquationAtPoint(const Point & p) const override;

  /**
   * @brief Return the polygon attributes for this object.
   *
   * @return map containing: num_sides (int), apothem (Real)
   */
  std::unordered_map<std::string, AttributeVariant> getAttributes() const override;

  // helper getter functions to get additional polygon dimensions/attributes

  /**
   * @brief Get the number of sides for the polygon
   *
   * @return number of sides
   */
  int getNumSides() { return _n_sides; }

  /**
   * @brief Get the polygon apothem (center-to-flat distance)
   *
   * @return apothem value
   */
  Real getApothem() { return _apothem; }

  /**
   * @brief Get the polygon side length
   *
   * @return side length
   */
  Real getSideLength() { return 2.0 * _apothem * std::tan(M_PI / _n_sides); }

  /**
   * @brief Get the the circumradius for the polygon (center-to-vertex distance)
   *
   * @return radius value
   */
  Real getRadius() { return _apothem / (std::cos(M_PI / _n_sides)); }

protected:
  /**
   * @brief Return a deep copy of this unit.
   *
   * @return unique_ptr to a new CSGNPolygonUnit with identical parameters
   */
  std::unique_ptr<CSGSurface> clone() const override
  {
    return std::make_unique<CSGNPolygonUnit>(_name, _n_sides, _apothem);
  }

  /**
   * @brief Create N CSGPlane surfaces in _internal_base, one per polygon side.
   *
   * Stores const pointers to the registered surfaces for use in getExpandedRegion().
   */
  void expandUnit() override;

private:
  /// Number of sides of the regular polygon
  const int _n_sides;

  /// Distance from the polygon center to the midpoint of each side
  const Real _apothem;

  /// Pointers to the expanded plane surfaces, populated during expandUnit()
  std::vector<const CSGSurface *> _planes;

#ifdef MOOSE_UNIT_TEST
  FRIEND_TEST(CSGEngUnitTest, testPolygonUnitExpansion);
  FRIEND_TEST(CSGEngUnitTest, testPolygonUnitClone);
#endif
};

} // namespace CSG

#include "CSGNPolygonUnit.h"
#include "CSGPlane.h"
#include "CSGBase.h" // complete type required to call methods on _internal_base

namespace CSG
{

CSGNPolygonUnit::CSGNPolygonUnit(const std::string & name, int n_sides, Real apothem)
  : CSGSurfaceEngUnit(name), _n_sides(n_sides), _apothem(apothem)
{
  if (_n_sides < 3)
    mooseError("N-sided polygon engineering unit " + name + " must have 3 or more sides.");
  if (_apothem <= 0.0)
    mooseError("N-sided polygon engineering unit " + name + " apothem must be positive.");
}

Real
CSGNPolygonUnit::evaluateSurfaceEquationAtPoint(const Point & p) const
{
  // Condition for interior (negative value) vs exterior (positive value) for the polygon is:
  // for point (x, y), if the following is true for all values of k, then the point is "interior".
  // If any one value is positive, then the point is outside the polygon. Therefore, we calculate
  // the maximum value. If that max value remains negative, then getHalfspaceFromPoint will
  // correctly determine the point to be interior.
  //
  //    x*cos(2*pi*k/N) + y*sin(2*pi*k/N) <= apothem, for all k = 0..N-1

  Real max_val =
      -std::numeric_limits<Real>::max(); // initialize to extremely large negative (to be updated)
  for (int k = 0; k < _n_sides; ++k)
  {
    auto val =
        (p(0) * std::cos(2.0 * M_PI * k / _n_sides) + p(1) * std::sin(2.0 * M_PI * k / _n_sides)) -
        _apothem;
    if (val > max_val)
      max_val = val;
  }
  return max_val;
}

std::unordered_map<std::string, AttributeVariant>
CSGNPolygonUnit::getAttributes() const
{
  return {{"num_sides", _n_sides}, {"apothem", _apothem}};
}

void
CSGNPolygonUnit::expandUnit()
{
  // Polygon orientation assumes an infinite prism oriented with the z-axis.
  // The right-most face is parallel to the y-axis and is centered at the origin.
  // Equation for the kth face, where the 0th face is the right-most and A is the apothem, follows
  // this equation:
  //    x*cos(2*pi*k/N) + y*sin(2*pi*k/N) + z*0 = A
  // Coefficients for the plane ax + by + cz = d:
  //    a = cos(2*pi*k/N)
  //    b = sin(2*pi*k/N)
  //    c = 0.0
  //    d = A (apothem)
  //
  // Surface naming scheme: [UnitName]_expanded_surf_[k]

  Real a, b; // to be calculated based on side
  Real c = 0.0;
  Real d = _apothem;

  // Initialize region to be added to:
  Point p(0, 0, 0); // origin used for determining half-space

  // base name for surfaces
  std::string base_name = getName() + "_expanded_surf_";

  for (int k = 0; k < _n_sides; ++k)
  {
    auto sname = base_name + std::to_string(k);
    a = std::cos(2.0 * M_PI * k / _n_sides);
    b = std::sin(2.0 * M_PI * k / _n_sides);
    std::unique_ptr<CSG::CSGPlane> s_ptr = std::make_unique<CSG::CSGPlane>(sname, a, b, c, d);
    auto & surf = _internal_base->addSurface(std::move(s_ptr));

    // determine the half-space that contains the origin for this surface
    auto hp_type = surf.getHalfspaceFromPoint(p);
    // half-space region for this surface only (to be intersected below)
    auto hp = (hp_type == CSGSurface::Halfspace::POSITIVE) ? +surf : -surf;

    // start the region with first half-space, otherwise intersect with existing region
    if (_expanded_region.getRegionType() == CSGRegion::RegionType::EMPTY)
      _expanded_region = hp;
    else
      _expanded_region &= hp; // intersect with existing region
  }
}

} // namespace CSG