Implicit Complement and “Graveyard” in OpenMC

With the constants addressed, the Implicit Complement works well in OpenMC. I settled the descrepancy after adding the graveyard as well. All material IDs were being set (incorrectly) on the Fortran side of the code.

I also updated the location for the next_cell call to DAGMC to come after OpenMC makes many decisions about how the particle should behave when crossing a surface. It is called only once the next volume is searched for. In the future, I’d like to hack surfaces to have positive and negative neighbor lists in OpenMC which are just their volumes with a forward and reverse sense, respectively.

After adding more boundary condition options for surfaces and forwarding more geometry calls to DAGMC, this is pretty much ready to go!

PyMOAB in PyNE

I started looking back at this work again as well. The place I left off was when I discovered that MOAB’s ScdInterface has no robust way to recover structured mesh sets when loading a file. These sets are tagged with a BOX_DIMS tag, but other than that no information about the ordering, starting point, etc. of the mesh is stored.

So how was PyNE doing this?

PyNE relies on elements to come back in an ijk ordering when using a PyTAPs MeshSet method called iterate. My instinct tells me this eventually backs onto a MOAB interface call named create_set_iterator (or the like). Here is the function signature for that call in MOAB:

    /** \brief Create an iterator over the set
     * Create a new iterator that iterates over entities with the specified type or dimension.
     * Only one of ent_type or dim can be set; use dim=-1 or ent_type=MBMAXTYPE for the other.
     * Iterators for list-type (ordered) sets are stable over set modification, unless entity
     * removed or deleted is the one at the current position of the iterator.  If the check_valid
     * parameter is passed as true, entities are checked for validity before being passed back by
     * get_next_entities function (checking entity validity can have a non-negligible cost).
     *
     * Iterators returned by this function can be deleted using the normal C++ delete function.
     * After creating the iterator through this function, further interactions are through methods
     * on the SetIterator class.
     * \param meshset The entity set associated with this iterator (use 0 for whole instance)
     * \param ent_type Entity type associated with this iterator
     * \param ent_dim Dimension associated with this iterator
     * \param chunk_size Chunk size of the iterator
     * \param check_valid If true, entities are checked for validity before being returned
     */
  virtual ErrorCode create_set_iterator(EntityHandle meshset,
                                        EntityType ent_type,
                                        int ent_dim,
                                        int chunk_size,
                                        bool check_valid,
                                        SetIterator *&set_iter);

I need to trace this back though to find out how PyTAPS/iTAPs actually makes this call and forwards information. I’ll then need to add this capability to PyMOAB in order to be able to continue work on this.

From what I can see in the PyTAPS source code, the MeshSet.iterate method definition is in iMesh_entSet.inl in PyTAPS. When the iterate method is called, this forwards to the iMesh interface PyTAPS has been built upon. In our case this is MOAB’s iMesh interface. Here, there are two iterator classes, iBase_EntArrIterator and iBase_EntIterator, which are used to iterate over sets of entities in the underlying mesh system the iMesh interface is linked too. It turns out that when these are created, it does not call the MOAB core function I expected it to. Instead, a set of rules are followed for the set of entities returned. These entities are then stored in a container (either a std::vector or a moab::Range) and returned as need through the iBase iterator classes I’ve mentioned above. The source for one of the iterator classes is listed here:

template <class Container> class MBIter : public iBase_EntityArrIterator_Private
{
  protected:
    Container iterData;
    typename Container::const_iterator iterPos;

  public:
    MBIter( iBase_EntityType type,
            iMesh_EntityTopology topology,
            EntityHandle set,
            int arr_size,
            bool recursive = false)
            : iBase_EntityArrIterator_Private( type, topology, set, arr_size, recursive ),
        iterPos(iterData.end()) {}

    ~MBIter() {}

    typename Container::const_iterator position() const {return iterPos;};

    typename Container::const_iterator end() const {return iterData.end();};

    ErrorCode step(int num_steps, bool &at_end)
    {
      return step_iterator(iterPos, end(), num_steps, at_end);
    }

    void get_entities( Core* mb, EntityHandle* array, int& count )
    {
      for (count = 0; count < arrSize && iterPos != iterData.end(); ++iterPos)
        if (mb->is_valid(*iterPos))
          array[count++] = *iterPos;
    }

    virtual ErrorCode reset( Interface* mb ) {
      ErrorCode result;
      iterData.clear();
      if (entTopo != iMesh_ALL_TOPOLOGIES) {
        if (entTopo == iMesh_SEPTAHEDRON)
          result = MB_SUCCESS;
        else
          result = mb->get_entities_by_type( entSet, mb_topology_table[entTopo], iterData, isRecursive );
      }
      else if (entType != iBase_ALL_TYPES) {
        result = mb->get_entities_by_dimension( entSet, entType, iterData, isRecursive );
        if (entType == iBase_REGION)
          remove_type( iterData, MBKNIFE );
      }
      else {
        result = mb->get_entities_by_handle( entSet, iterData, isRecursive );
        remove_type( iterData, MBENTITYSET );
        remove_type( iterData, MBKNIFE );
      }
      iterPos = iterData.begin();
      return result;
    }
};