import from dgcode git: version 3.3.1

CHANGELOG

@@ -1,3 +1,17 @@
+v3.3.1 2022-10-06
+      * Fix info.getComposition() for multiphase materials in the Python
+        interface.  This issue is not believed to have affected many users, but
+        a fix is crucial in order to support multiphase materials in OpenMC
+        hooks which are currently being prepared (github issue #102).
+      * Minor fix to make Geant4 hooks work with Geant4 release 11. Note that
+        this is still not the long awaited major update which makes the hooks
+        work with Geant4 in multithreaded mode.
+      * Rewrote internal code in NCIter.hh to avoid usage of deprecated
+        std::iterator. This avoids compilation warnings on newer platforms
+        (github issue #100).
+      * Reduce precision from %.15g to %.14g in output of ncrystal_ncmat2hkl to
+        get reproducible results across platforms.
+
 v3.3.0 2022-08-23
       * Updates to CMake code which is needed to support creation of
         NCrystal conda packages.

CMakeLists.txt

@@ -69,7 +69,7 @@ if( NOT NCRYSTAL_NOTOUCH_CMAKE_BUILD_TYPE )
 endif()
 
 #Setup project:
-project(NCrystal VERSION 3.3.0 LANGUAGES CXX C)
+project(NCrystal VERSION 3.3.1 LANGUAGES CXX C)
 
 if( NOT NCRYSTAL_NOTOUCH_CMAKE_BUILD_TYPE )
   if ("${CMAKE_BUILD_TYPE}" STREQUAL "")

VERSION

@@ -1 +1 @@
-3.3.0
+3.3.1

cmake/ncrystal-config.in

@@ -70,7 +70,7 @@ if installation_info.get('libdir','@').startswith('@'):
 def hasopt(optname):
     return installation_info['opt_%s'%optname]
 
-####Remove some directories if nothing was installed there due configuration options:
+####Remove some directories if nothing was installed there due to configuration options:
 ###if not ( hasopt('data') and not hasopt('embed_data') ):
 ###    installation_info['datadir']=None
 ###if not ( hasopt('mcstas') and not hasopt('embed_data') ):

docs/ncmat_doc.md

@@ -722,7 +722,7 @@ indicate the state of matter if needed.
 
 ## Changes for the @DYNINFO section ##
 
-In files without a @DEBYETEMPERATURE section, it is now allowd to specify the
+In files without a @DEBYETEMPERATURE section, it is now allowed to specify the
 Debye temperatures directly inside @DYNINFO sections of type *vdosdebye*,
 providing them in units of kelvin via the new keyword `debye_temp`, as for
 instance:

ncrystal_core/include/NCrystal/NCInfoTypes.hh

@@ -101,6 +101,9 @@ namespace NCrystal {
     static constexpr const char * unit() noexcept { return ""; }
     static AtomIndex createInvalidObject();
     constexpr bool isInvalid() const noexcept;
+    //NB: Since unsigned is only guaranteed to be at least 16bit, the underlying
+    //value representing invalid indices might be as low as 65535 (of course, in
+    //practice it is usually 4294967295).
   };
 
   struct NCRYSTAL_API IndexedAtomData {

ncrystal_core/include/NCrystal/NCVersion.hh

@@ -39,9 +39,9 @@
 
 #define NCRYSTAL_VERSION_MAJOR 3
 #define NCRYSTAL_VERSION_MINOR 3
-#define NCRYSTAL_VERSION_PATCH 0
-#define NCRYSTAL_VERSION   3003000 /* (1000000*MAJOR+1000*MINOR+PATCH)   */
-#define NCRYSTAL_VERSION_STR "3.3.0"
+#define NCRYSTAL_VERSION_PATCH 1
+#define NCRYSTAL_VERSION   3003001 /* (1000000*MAJOR+1000*MINOR+PATCH)   */
+#define NCRYSTAL_VERSION_STR "3.3.1"
 
 #include "NCrystal/ncapi.h"
 #include <stdexcept>

ncrystal_core/include/NCrystal/internal/NCIter.hh

@@ -22,8 +22,6 @@
 ////////////////////////////////////////////////////////////////////////////////
 
 #include "NCrystal/NCDefs.hh"
-#include <iterator>
-#include <initializer_list>
 
 namespace NCrystal {
 
@@ -33,51 +31,58 @@ namespace NCrystal {
   //
   //Examples:
   //
-  //for (auto&& e : enumerate(mycontainer))
+  //for ( auto&& e : enumerate(mycontainer) )
   //   othercontainer[e.idx] = e.val;
-  //for (auto&& e : enumerate(mycontainer))
+  //for ( auto&& e : enumerate(mycontainer) )
+  //   e.val += 0.1*e.idx;
+  //
+  //The by-value form is also OK, but it might still be possible to modify the
+  //elements:
+  //
+  //for ( auto e : enumerate(mycontainer) )
   //   e.val += 0.1*e.idx;
   //
   //It is also possible to explicitly impose read-only access (for code clarity perhaps):
   //
   //  for (const auto& e : enumerate(mycontainer))
   //     othercontainer[e.idx] = e.val;
   //
-  //The form (auto& e: enumerate(mycontainer)) does not always compile
-  //(depending on constness), but is anyway not needed given the options above.
+  // The form (auto& e: enumerate(mycontainer)) does not always compile
+  // (depending on constness), but is anyway not needed given the options above.
+
+  namespace detail {
+    template <typename TIter>
+    struct EnumIter
+    {
+      std::size_t idx;
+      TIter it;
+      bool operator != (const EnumIter & other) const { return it != other.it; }
+      void operator ++() { ++idx; ++it; }
+      struct Entry {
+        std::size_t idx;
+        decltype( *TIter()) val;
+      };
+      Entry operator*() const { return Entry{idx, *it}; }
+    };
+
+    template <typename T, typename TIter>
+    struct EnumIterWrapper
+    {
+      static_assert(std::is_reference<T>::value,"Inefficient enumerate(..) usage");
+      T container;
+      EnumIter<TIter> begin() { return EnumIter<TIter>{ 0, std::begin(container) }; }
+      EnumIter<TIter> end() { return EnumIter<TIter>{ 0, std::end(container) }; }
+    };
 
-  template <typename TIter>
-  class enumerate_iter_t : std::iterator<std::forward_iterator_tag, typename std::iterator_traits<TIter>::value_type> {
-  public:
-    using reference = typename std::iterator_traits<TIter>::reference;
-  private:
-    std::size_t m_idx = 0;
-    TIter m_it;
-  public:
-    explicit enumerate_iter_t(TIter&& it): m_it{ std::forward<TIter>(it) } {}
-    enumerate_iter_t& operator++() { ++m_it; ++m_idx; return *this; }
-    bool operator!=(const enumerate_iter_t& o) const { return m_it != o.m_it; }
-    struct entry_t { const std::size_t idx; reference val; };
-    struct const_entry_t { const std::size_t idx; const reference val; };
-    entry_t operator*() { return { m_idx, *m_it }; }
-    const_entry_t operator*() const { return { m_idx, *m_it }; }
-  };
-  template <typename T> struct enumerate_wrapper_t { T& container; };
-  template <typename T>
-  auto begin(enumerate_wrapper_t<T>& w) -> enumerate_iter_t<decltype(std::begin(w.container))>
-  {
-    return enumerate_iter_t<decltype(std::begin(w.container))>(std::begin(w.container));
   }
-  template <typename T>
-  auto end(enumerate_wrapper_t<T>& w) -> enumerate_iter_t<decltype(std::begin(w.container))>
+  template <typename T,
+            typename TIter = decltype(std::begin(std::declval<T>())),
+            typename TIterE = decltype(std::end(std::declval<T>()))>
+  inline constexpr detail::EnumIterWrapper<T,TIter> enumerate( T && container )
   {
-    return enumerate_iter_t<decltype(std::begin(w.container))>(std::end(w.container));
+    static_assert(std::is_same<TIter,TIterE>::value,"");
+    return detail::EnumIterWrapper<T,TIter>{ std::forward<T>(container) };
   }
-  template <typename T>
-  enumerate_wrapper_t<T> enumerate(T&& container) { return {container}; }
-  template <typename T>
-  enumerate_wrapper_t<std::initializer_list<T>> enumerate(std::initializer_list<T>&& container) { return {container}; }
-
 }
 
 

ncrystal_core/include/NCrystal/ncrystal.h

@@ -274,7 +274,10 @@ extern "C" {
   NCRYSTAL_API double ncrystal_info_dspacing_from_hkl( ncrystal_info_t, int h, int k, int l );
 
 
-  /* Composition (always >=1 component):                                           */
+  /* Composition (always >=1 component). Note that for multiphase objects, the     */
+  /* provided atomdataidx will be invalid, so it is important in general to use    */
+  /* ncrystal_create_component_atomdata(..) to access the atomdata object          */
+  /* associated with the i'th component, and NOT ncrystal_create_atomdata(..):     */
   NCRYSTAL_API unsigned ncrystal_info_ncomponents( ncrystal_info_t );
   NCRYSTAL_API void ncrystal_info_getcomponent( ncrystal_info_t, unsigned icomponent,
                                                 unsigned* atomdataindex,
@@ -288,6 +291,10 @@ extern "C" {
   NCRYSTAL_API ncrystal_atomdata_t ncrystal_create_atomdata( ncrystal_info_t,
                                                              unsigned atomdataindex );
 
+  /* Same, but via index in composition vector:                                    */
+  NCRYSTAL_API ncrystal_atomdata_t ncrystal_create_component_atomdata( ncrystal_info_t,
+                                                                       unsigned icomponent );
+
   /* Get atom data fields. Each object falls in one of three categories:           */
   /* 1) Natural elements (ncomponents=A=0,Z>0)                                     */
   /* 2) Single isotope (ncomponents=0, Z>0, A>=Z)                                  */
@@ -516,9 +523,9 @@ extern "C" {
 #endif
 #define NCRYSTAL_VERSION_MAJOR 3
 #define NCRYSTAL_VERSION_MINOR 3
-#define NCRYSTAL_VERSION_PATCH 0
-#define NCRYSTAL_VERSION   3003000 /* (1000000*MAJOR+1000*MINOR+PATCH)             */
-#define NCRYSTAL_VERSION_STR "3.3.0"
+#define NCRYSTAL_VERSION_PATCH 1
+#define NCRYSTAL_VERSION   3003001 /* (1000000*MAJOR+1000*MINOR+PATCH)             */
+#define NCRYSTAL_VERSION_STR "3.3.1"
   NCRYSTAL_API int ncrystal_version(); /* returns NCRYSTAL_VERSION                  */
   NCRYSTAL_API const char * ncrystal_version_str(); /* returns NCRYSTAL_VERSION_STR */
 

ncrystal_core/src/ncrystal.cc

@@ -50,11 +50,12 @@ namespace NCrystal {
     class AtomDataObj : private MoveOnly {
       //Hold AtomData and related info. This wrapper class is associated with a
       //particular Info instance, and therefore also knows the displayLabel
-      //(unless a sub-component). It also caches some strings which would
-      //otherwise be temporary objects, making it safe to return their values to
-      //C code. Note that the displayLabels are only valid in connection with a
-      //particular Info object, and only "top-level" atoms, the ones with an
-      //AtomIndex on the Info object, have a displayLabel.
+      //(unless a sub-component, or atomdata is from the composition vector). It
+      //also caches some strings which would otherwise be temporary objects,
+      //making it safe to return their values to C code. Note that the
+      //displayLabels are only valid in connection with a particular Info
+      //object, and only "top-level" atoms, the ones with an AtomIndex on the
+      //Info object, have a displayLabel.
       shared_obj<const AtomData> m_atomDataSO;
       std::unique_ptr<std::string> m_displayLabel_ptr;
       std::unique_ptr<std::string> m_description_ptr;
@@ -1439,8 +1440,26 @@ ncrystal_atomdata_t ncrystal_create_atomdata( ncrystal_info_t ci,
 {
   try {
     auto& info = ncc::extract(ci);
-    return ncc::createNewCHandle<ncc::Wrapped_AtomData>( info->atomDataSP(NC::AtomIndex{atomdataindex}),
-                                                         info->displayLabel(NC::AtomIndex{atomdataindex}) );
+    NC::AtomIndex aidx{atomdataindex};
+    if ( aidx.isInvalid() )
+      NCRYSTAL_THROW2(BadInput,"ncrystal_create_atomdata: provided atomdataidx is invalid.");
+    return ncc::createNewCHandle<ncc::Wrapped_AtomData>( info->atomDataSP(aidx),
+                                                         info->displayLabel(aidx) );
+  } NCCATCH;
+  return {nullptr};
+}
+
+ncrystal_atomdata_t ncrystal_create_component_atomdata( ncrystal_info_t nfo,
+                                                        unsigned icomponent )
+{
+  try {
+    auto& info = ncc::extract(nfo);
+    auto n = info->getComposition().size();
+    nc_assert(n>=1);
+    if ( ! ( icomponent<n ) )
+      NCRYSTAL_THROW(BadInput,"Requested component index is out of bounds");
+    const auto& comp = info->getComposition().at(icomponent);
+    return ncc::createNewCHandle<ncc::Wrapped_AtomData>( comp.atom.atomDataSP );
   } NCCATCH;
   return {nullptr};
 }
@@ -1501,9 +1520,10 @@ void ncrystal_info_getcomponent( ncrystal_info_t ci, unsigned icomponent,
     auto& info = ncc::extract(ci);
     auto n = info->getComposition().size();
     nc_assert(n>=1);
-    if ( ! ( icomponent<n) )
+    if ( ! ( icomponent<n ) )
       NCRYSTAL_THROW(BadInput,"Requested component index is out of bounds");
     const auto& comp = info->getComposition().at(icomponent);
+    //NB: comp will contain invalid atom index on multiphase object!
     *atomdataindex = comp.atom.index.get();
     *fraction = comp.fraction;
     return;

ncrystal_geant4/include/G4NCrystal/G4NCManager.hh

@@ -131,8 +131,14 @@ namespace G4NCrystal {
     //we can't be sure it is always the same in different materials. Unfortunately
     //GetConstProperty triggers a G4 exception in case the key is not present,
     //so for safety we have to trigger one extra map search via ConstPropertyExists:
+
+#  if G4VERSION_NUMBER >= 1100
+    if ( matprop->GetMaterialConstPropertyNames().empty() || !matprop->ConstPropertyExists(m_key) )
+      return not_found;
+#  else
     if ( matprop->GetConstPropertyMap()->empty() || !matprop->ConstPropertyExists(m_key) )
       return not_found;
+#endif
     return static_cast<unsigned>(matprop->GetConstProperty(m_key));
 #endif
   }

ncrystal_geant4/src/G4NCManager.cc

@@ -105,7 +105,11 @@ void NCG4::Manager::addScatterProperty(G4Material* mat,NCrystal::ProcImpl::ProcP
     idx = it->second;
   }
   assert( unsigned(double(idx)) == idx );//make sure we can get the idx back out
-  matprop->AddConstProperty(m_key.c_str(), idx);
+  matprop->AddConstProperty(m_key.c_str(), idx
+#if G4VERSION_NUMBER >= 1100
+                            , true//createNewKey=true is required explicitly since G4 11.0
+#endif
+                            );
 }
 
 void NCG4::Manager::cleanup()

ncrystal_mcstas/NCrystal_sample.comp

@@ -21,7 +21,7 @@
 *
 * %I
 * Written by: NCrystal developers
-* Version: 3.3.0
+* Version: 3.3.1
 * Origin: NCrystal Developers (European Spallation Source ERIC and DTU Nutech)
 *
 * McStas sample component for the NCrystal library for thermal neutron transport

ncrystal_python/__init__.py

@@ -52,7 +52,7 @@
 
 #NB: Synchronize meta-data below with fields in setup.py.in meta data:
 __license__ = "Apache 2.0, http://www.apache.org/licenses/LICENSE-2.0"
-__version__ = '3.3.0'
+__version__ = '3.3.1'
 __status__ = "Production"
 __author__ = "NCrystal developers (Thomas Kittelmann, Xiao Xiao Cai)"
 __copyright__ = "Copyright 2015-2022 %s"%__author__
@@ -464,6 +464,8 @@ def ncrystal_info_getcomp(nfo,icomp):
     _raw_atomdata_getfields=_wrap('ncrystal_atomdata_getfields',None,(ncrystal_atomdata_t,_cstrp,_cstrp,
                                                                       _dblp,_dblp,_dblp,_dblp,
                                                                       _uintp,_uintp,_uintp),hide=True)
+    _wrap('ncrystal_create_component_atomdata',ncrystal_atomdata_t,(ncrystal_info_t,_uint))
+
     def ncrystal_atomdata_createsubcomp(ad,icomp):
         fraction = _dbl()
         comp_ad = _raw_atomdata_subcomp(ad,icomp,fraction)
@@ -977,8 +979,11 @@ def displayLabel(self):
         return self.__dl
 
     def isTopLevel(self):
-        """Whether or not AtomData appears directly on an Info object (if not, it must
-        be a component (direct or indirect) of a top level AtomData object"""
+        """Whether or not AtomData appears directly on an Info object. If not,
+         it will most likely either be a component (direct or indirect) of a top
+         level AtomData object, or be taken from the composition list of a
+         multi-phase object.
+        """
         return bool(self.__dl)
 
     def description(self,includeValues=True):
@@ -1084,10 +1089,22 @@ def getPhases(self):
     def _initComp(self):
         assert self.__comp is None
         nc = _rawfct['ncrystal_info_ncomponents'](self._rawobj)
-        self.__comp = []
+        l = []
         for icomp in range(nc):
             atomidx,fraction = _rawfct['ncrystal_info_getcomp'](self._rawobj,icomp)
-            self.__comp += [(fraction,self._provideAtomData(atomidx))]
+            #NB: atomidx will be invalid in case of multiphase objects!
+            if atomidx < 65535:
+                #Most likely a single-phase object with valid atomidx, we can
+                #use self._provideAtomData and share the AtomData objects also here on the python side:
+                l += [(fraction,self._provideAtomData(atomidx))]
+            else:
+                #Most likely a multi-phase object with invalid atomidx, we must
+                #create new AtomData objects, based on ncrystal_create_component_atomdata:
+                raw_ad = _rawfct['ncrystal_create_component_atomdata'](self._rawobj,icomp)
+                obj = AtomData(raw_ad)
+                assert not obj.isTopLevel()#does not appear directly on Info object
+                l += [(fraction,obj)]
+        self.__comp = l
         return self.__comp
 
     def stateOfMatter(self):
@@ -1238,7 +1255,8 @@ def getStructureInfo(self):
 
     def _provideAtomData(self,atomindex):
         if atomindex >= len(self.__atomdatas):
-            assert atomindex < 100000#sanity check
+            if atomindex >= 65535:
+                raise NCLogicError(f'Invalid atomindex ({atomindex}) provided to Info._provideAtomData')
             self.__atomdatas.extend([None,]*(atomindex+1-len(self.__atomdatas)))
         obj = self.__atomdatas[atomindex]
         if obj:

ncrystal_python/mcstasutils.py

@@ -149,7 +149,7 @@ def cfgstr_2_union_instrument_code( *, cfgstr, name, split_by_physics = False ):
 
     return res
 
-def cfgstr_2_hkl(*, cfgstr, tgtformat, verbose=True, fp_format = '%.15g'):
+def cfgstr_2_hkl(*, cfgstr, tgtformat, verbose=True, fp_format = '%.14g'):
     """Function which can be used to create input files with reflections for
     McStas crystalline sample components like PowderN and Single_crystal, based
     on NCrystal cfg-strings (usually referring to NCMAT files with crystalline