From 77d99d499ba7a2b91e2145d9735f34c99351b61d Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 6 Aug 2020 14:19:53 -0500
Subject: [PATCH] Allow access to material temperature from openmc.lib

---
 include/openmc/material.h | 14 +++++++++-----
 openmc/lib/material.py    | 15 +++++++++++++++
 src/geometry_aux.cpp      | 12 ++----------
 src/material.cpp          | 18 ++++++++++++++++++
 4 files changed, 44 insertions(+), 15 deletions(-)

diff --git a/include/openmc/material.h b/include/openmc/material.h
index 93a90d9e4..9adb9ff97 100644
--- a/include/openmc/material.h
+++ b/include/openmc/material.h
@@ -131,6 +131,10 @@ public:
   //! \return Volume in [cm^3]
   double volume() const;
 
+  //! Get temperature of material
+  //! \return Temperature in [K]
+  double temperature() const;
+
   //----------------------------------------------------------------------------
   // Data
   int32_t id_ {C_NONE}; //!< Unique ID
@@ -154,11 +158,6 @@ public:
   // Thermal scattering tables
   std::vector<ThermalTable> thermal_tables_;
 
-  //! \brief Default temperature for cells containing this material.
-  //!
-  //! A negative value indicates no default temperature was specified.
-  double temperature_ {-1};
-
   std::unique_ptr<Bremsstrahlung> ttb_;
 
 private:
@@ -180,6 +179,11 @@ private:
   //----------------------------------------------------------------------------
   // Private data members
   gsl::index index_;
+
+  //! \brief Default temperature for cells containing this material.
+  //!
+  //! A negative value indicates no default temperature was specified.
+  double temperature_ {-1};
 };
 
 //==============================================================================
diff --git a/openmc/lib/material.py b/openmc/lib/material.py
index 0fbce7245..f3c794778 100644
--- a/openmc/lib/material.py
+++ b/openmc/lib/material.py
@@ -35,6 +35,9 @@ _dll.openmc_material_get_densities.errcheck = _error_handler
 _dll.openmc_material_get_density.argtypes = [c_int32, POINTER(c_double)]
 _dll.openmc_material_get_density.restype = c_int
 _dll.openmc_material_get_density.errcheck = _error_handler
+_dll.openmc_material_get_temperature.argtypes = [c_int32, POINTER(c_double)]
+_dll.openmc_material_get_temperature.restype = c_int
+_dll.openmc_material_get_temperature.errcheck = _error_handler
 _dll.openmc_material_get_volume.argtypes = [c_int32, POINTER(c_double)]
 _dll.openmc_material_get_volume.restype = c_int
 _dll.openmc_material_get_volume.errcheck = _error_handler
@@ -86,6 +89,12 @@ class Material(_FortranObjectWithID):
         List of nuclides in the material
     densities : numpy.ndarray
         Array of densities in atom/b-cm
+    name : str
+        Name of the material
+    temperature : float
+        Temperature of the material in [K]
+    volume : float
+        Volume of the material in [cm^3]
 
     """
     __instances = WeakValueDictionary()
@@ -141,6 +150,12 @@ class Material(_FortranObjectWithID):
         name_ptr = c_char_p(name.encode())
         _dll.openmc_material_set_name(self._index, name_ptr)
 
+    @property
+    def temperature(self):
+        temperature = c_double()
+        _dll.openmc_material_get_temperature(self._index, temperature)
+        return temperature.value
+
     @property
     def volume(self):
         volume = c_double()
diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp
index c58af934b..eff7016ab 100644
--- a/src/geometry_aux.cpp
+++ b/src/geometry_aux.cpp
@@ -181,17 +181,9 @@ assign_temperatures()
       if (i_mat == MATERIAL_VOID) {
         // Set void region to 0K.
         c->sqrtkT_.push_back(0);
-
       } else {
-        if (model::materials[i_mat]->temperature_ >= 0) {
-          // This material has a default temperature; use that value.
-          auto T = model::materials[i_mat]->temperature_;
-          c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * T));
-        } else {
-          // Use the global default temperature.
-          c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN *
-            settings::temperature_default));
-        }
+        const auto& mat {model::materials[i_mat]};
+        c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature()));
       }
     }
   }
diff --git a/src/material.cpp b/src/material.cpp
index a77d51ed8..1d4e18108 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -987,6 +987,12 @@ double Material::volume() const
   return volume_;
 }
 
+double Material::temperature() const
+{
+  // If material doesn't have an assigned temperature, use global default
+  return temperature_ >= 0 ? temperature_ : settings::temperature_default;
+}
+
 void Material::to_hdf5(hid_t group) const
 {
   hid_t material_group = create_group(group, "material " + std::to_string(id_));
@@ -1339,6 +1345,18 @@ openmc_material_get_id(int32_t index, int32_t* id)
   }
 }
 
+extern "C" int
+openmc_material_get_temperature(int32_t index, double* temperature)
+{
+  if (index < 0 || index >= model::materials.size()) {
+    set_errmsg("Index in materials array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  *temperature = model::materials[index]->temperature();
+  return 0;
+}
+
+
 extern "C" int
 openmc_material_get_volume(int32_t index, double* volume)
 {