simplify mechanism to detect if geometry entity is DAG

[gridley]

Description

There is a somewhat cumbersome CSGSurface class that does nothing aside from set a variable called geom_type. We have written accessors, stored this variable on a few types of geometry class, and defined an enum just to handle distinguishing if a surface is a DAGSurface.

Since this variable geom_type is only a GeometryType::DAG if a given entity is of just one subclass type (DAGSurface), it makes more sense to define a virtual method to handle this. Virtuals are used extensively elsewhere on these geometric entities, so this seamlessly solves the problem of checking whether a surface is a DAGSurface. This is also applied to the cell and universe classes.

I ended up making these changes because I want to play with defining a new type of surface, and the distinction of CSGSurface didn't make any sense to me here. In fact, I would argue all surfaces here are CSG surfaces; it would be more apt to call this class a NonDAGSurface. And of course, that would be silly. It's cleanest to just delete.

Fixes # (issue)

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 15) on any C++ source files (if applicable)
• [] I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[pshriwise]

I'm on board with most of the changes here @gridley!

I think you're right, the CSGSurface is redundant and causes the code to be overly verbose. The way that the surface types inherit from surface makes sense conceptually too -- a DAGMC surface of triangles is just another type of surface.

I would ask that we stick with the GeometryType enumerator though. There are types of geometry (e.g. volumetric mesh) that are on their way and the enumerator is more extensible than a specific method in the case that more geometry types/libraries crop up in the future.

[gridley]

Hm, gotcha. Can you give an example where building out conditionals on the different enum cases makes more sense than a virtual method? I would think that if we sufficiently design the base classes and behaviors, you shouldn't need to toggle different behaviors at the level of the calling code. For example maybe there should just be a cross_cell method on the universe, rather than externally switching between different cell search methods in geometry.cpp. I get the feeling that if we are having to do that kind of switching outside, the class could have been designed better and we are suffering from the Fortran hangover still.

But maybe I'm missing an example where this truly makes sense to switch among different behaviors at the outer calling code level.

Take this for example:

void Particle::cross_surface(const Surface& surf)
{

  if (settings::verbosity >= 10 || trace()) {
    write_message(1, "    Crossing surface {}", surf.id_);
  }

// if we're crossing a CSG surface, make sure the DAG history is reset
#ifdef DAGMC
  if (surf.geom_type_ == GeometryType::CSG)
    history().reset();
#endif

  // Handle any applicable boundary conditions.
  if (surf.bc_ && settings::run_mode != RunMode::PLOTTING) {
    surf.bc_->handle_particle(*this, surf);
    return;
  }

  // ==========================================================================
  // SEARCH NEIGHBOR LISTS FOR NEXT CELL

#ifdef DAGMC
  // in DAGMC, we know what the next cell should be
  if (surf.geom_type_ == GeometryType::DAG) {
    int32_t i_cell = next_cell(std::abs(surface()), cell_last(n_coord() - 1),
                       lowest_coord().universe) -
                     1;
    // save material and temp
    material_last() = material();
    sqrtkT_last() = sqrtkT();
    // set new cell value
    lowest_coord().cell = i_cell;
    auto& cell = model::cells[i_cell];

    cell_instance() = 0;
    if (cell->distribcell_index_ >= 0)
      cell_instance() = cell_instance_at_level(*this, n_coord() - 1);

    material() = cell->material(cell_instance());
    sqrtkT() = cell->sqrtkT(cell_instance());
    return;
  }
#endif

  bool verbose = settings::verbosity >= 10 || trace();
  if (neighbor_list_find_cell(*this, verbose)) {
    return;
  }

  // ==========================================================================
  // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS

  // Remove lower coordinate levels
  n_coord() = 1;
  bool found = exhaustive_find_cell(*this, verbose);

  if (settings::run_mode != RunMode::PLOTTING && (!found)) {
    // If a cell is still not found, there are two possible causes: 1) there is
    // a void in the model, and 2) the particle hit a surface at a tangent. If
    // the particle is really traveling tangent to a surface, if we move it
    // forward a tiny bit it should fix the problem.

    surface() = 0;
    n_coord() = 1;
    r() += TINY_BIT * u();

    // Couldn't find next cell anywhere! This probably means there is an actual
    // undefined region in the geometry.

    if (!exhaustive_find_cell(*this, verbose)) {
      mark_as_lost("After particle " + std::to_string(id()) +
                   " crossed surface " + std::to_string(surf.id_) +
                   " it could not be located in any cell and it did not leak.");
      return;
    }
  }
}

This is just not a great design. The universe itself should have a cross_surface method that acts on a GeometryState object. Putting the different logic types on the outside of where we've already added abstraction makes sense. Ideally the universe class should fully encapsulate all the auxiliary operations such as resetting the DAGMC particle history for example.

I asked Claude to put a name to the concerns I have here in terms of standard OOP jargon. Here is its response, which I think is pretty correct.

Based on your comment, you're describing several key software design principles and architectural concerns:

Primarily, you're advocating for the Encapsulation principle and specifically pointing out a violation of Tell, Don't Ask. Rather than asking about the type of geometry and then executing different logic branches, the behavior should be encapsulated within the appropriate object.
The code also appears to violate the Open/Closed Principle (from SOLID) - the conditional branches suggest the code needs modification every time a new geometry type is added, rather than being open for extension through polymorphism.
You're suggesting a move toward Polymorphism over type-based conditionals, which aligns with the general OOP principle of favoring dynamic dispatch over explicit type checking.
The specific suggestion about moving the cross_surface method to the Universe class and having it operate on a GeometryState reflects the Single Responsibility Principle - the Universe class should encapsulate all geometry-related behaviors.

Your comment about "Fortran hangover" particularly resonates - it's pointing out that the procedural-style switching logic is at odds with object-oriented principles, suggesting the code hasn't fully embraced modern OOP design patterns.
The specific pattern you're advocating for sounds like it would better align with the Strategy Pattern, where different geometry types could implement their own crossing behavior, rather than having it handled through conditional logic in the client code.

[gridley]

And let me add: I'd be happy to include a refactor like that as part of this PR. Ideally, we shouldn't have to check is_DAG anywhere, based on reflecting on this a little more.

[pshriwise]

Hm, gotcha. Can you give an example where building out conditionals on the different enum cases makes more sense than a virtual method?

Sorry, I should've been more specific here. What I was asking is that the virtual method return a GeometryType enumerator instead of a Boolean value. The data members are set immediately based on the object type, so it's a little silly to have a data member at all. Something like the following should do I think:

class Surface {

  virtual GeometryType geom_type() const { return GeometryType::CSG; };

};

To answer your question more directly though, no, I don't think building out conditionals makes more sense than a virtual method (assuming that you mean a single virtual method to handle surface crossings on the Universe class), but that's not what is being implemented in the current change set and didn't seem to be proposed in the PR description. The update to an is_DAG method still requires a branch in Particle::cross_surface and conditionals/methods would need to be built out for different geometry types as well -- we'd need to implement new methods (is_volumetric_mesh for example) to support new geometry types. So, at that level, I don't see a difference in design from what exists now (based on the current changes and PR description).

The implementation of is_DAG would arguably allow us to cleanly encapsulate some more complicated logic behind checking whether or not the object is DAGMC geometry or not, but this isn't really the case -- demonstrated by the fact that the returns are hardcoded. In contrast, I find the return of a GeoemtryType to be more clear when performing checks and require fewer methods to search through and implement (assuming extension to more Geometry types). It also seems (to me) to make sense for a CSG surface to have a geometry_type method but less so an is_DAG method. No need to allude to DAGMC at all on an XPlane object I'd say.

Happy to discuss further and let me know if I missed the mark on addressing your points here.

This is just not a great design. The universe itself should have a cross_surface method that acts on a GeometryState object. Putting the different logic types on the outside of where we've already added abstraction makes sense. Ideally the universe class should fully encapsulate all the auxiliary operations such as resetting the DAGMC particle history for example.

There are a number of areas where design could be improved -- 100% agreed!! Such changes would move the PR in a rather new direction and so I'd prefer such a refactor be performed separately. It's a PR that I'd certainly support and be more than happy to review. In light of the current proposed changes focused on simplifying the inheritance and in the absence of such a refactor, I'd still prefer we keep the GeometryType enum for now if you're open to that.

Thanks to Claude.ai for the input! I'll be making good use of the term "Fortran hangover" often I think.

[gridley]

OK, this all makes sense then! Baby steps. I'll get your suggested way of doing it in pretty soon.

[gridley]

Hey @pshriwise, sorry to take forever to get your suggestion in. I believe the PR as it stands is fully in line with your suggestion.

[gridley]

Looks like the CI fail is maybe due to some ncrystal changes that came through recently?

+ /home/runner/ncrystal_inst/bin/ncrystal-config --setup
ncrystal-config: error: Missing or invalid arguments. Run with --help for usage instructions.
Error: Process completed with exit code 1.

[pshriwise]

Looks like the CI fail is maybe due to some ncrystal changes that came through recently?

Hmmm I'll check that out. Is this perhaps related to #3274 @tkittel?

Otherwise the changes look good to me, thanks @gridley!

[tkittel]

Absolutely right @pshriwise, it is related to #3274 ! I am not sure why this did not cause the CI tests in #3274 fail, but from looking at gha-install-ncrystal.sh I can clearly see the reason: It uses the ncrystal develop branch (updated this weekend with a release candidate for NCrystal 4.0.0) where the ncrystal-config --setup flag has been removed.

I will try to add a fix for this as well to #3274 (i.e. rip out gha-install-ncrystal.sh and simply pip install ncrystal instead, which is now a thing).

[pshriwise]

Thanks for the quick response @tkittel! If there's anything I can do to help with #3274 let me know, but please don't interrupt your vacation significantly.