Merge pull request #88 from sandialabs/feature/read_properties_from_exo

Feature/read properties from exodus

optimism/ReadExodusMesh.py

@@ -1,8 +1,8 @@
-import netCDF4
-
 from optimism.JaxConfig import *
 from optimism import Mesh
 from optimism import Interpolants
+import netCDF4
+import numpy as onp
 
 exodusToNativeTri6NodeOrder = np.array([0, 3, 1, 5, 4, 2])
 
@@ -28,7 +28,20 @@ def read_exodus_mesh(fileName):
         return Mesh.Mesh(coords, conns, simplexNodesOrdinals, basis, basis1d,
                          blocks, nodeSets, sideSets)
 
-
+
+def read_exodus_mesh_element_properties(fileName, varNames, blockNum=1):
+    varValues = []
+    with netCDF4.Dataset(fileName) as exData:
+        blockNames = _read_names_list(exData, 'eb_names')
+        if len(blockNames) > 1:
+            raise ValueError('Only single blocks are supported currently')
+
+        for name in varNames:
+            varValues.append(_read_block_variable_values(exData, blockNum - 1, name))
+
+    return np.vstack(varValues).T
+
+
 def _read_coordinates(exodusDataset):
     nNodes = len(exodusDataset.dimensions['num_nodes'])
     nDims = len(exodusDataset.dimensions['num_dim'])
@@ -77,6 +90,35 @@ def _read_blocks(exodusDataset):
     return conns, blocks
 
 
+def _read_block_variable_values(exodusDataset, blockOrdinal, variableName):
+    # read element variables currently on mesh
+    record = exodusDataset.variables['name_elem_var']
+    record.set_auto_mask(False)
+
+    propNamesInMesh = []
+    for n in range(record.shape[0]):
+        propNamesInMesh.append(''.join(onp.char.decode(record[n, :])))
+
+    # make sure the requested variable is there
+    try:
+        idx = propNamesInMesh.index(variableName)
+    except:
+        string = f'Requested variable {variableName} not found on mesh.\n'
+        string += 'Available variables are:\n'
+        for name in propNamesInMesh:
+            string += f'  {name}\n'
+        print(string)
+        raise KeyError(string)
+
+    # get the name correct
+    key = f'vals_elem_var{idx + 1}eb{blockOrdinal + 1}'
+
+    # finally read the variable
+    record = exodusDataset.variables[key]
+    record.set_auto_mask(False)
+    return np.array(record[:])
+
+
 def _read_node_sets(exodusDataset):
     if "num_node_sets" in exodusDataset.dimensions:
         nodeSetNames = _read_names_list(exodusDataset, "ns_names")

test/test_ReadExodusMesh.py

@@ -1,4 +1,5 @@
 import pathlib # to reach sample mesh file
+import pytest
 import sys
 import unittest
 
@@ -150,5 +151,25 @@ def objective(Uu):
         self.assertNear(np.linalg.norm(grad_func(Uu)), 0.0, 14)
 
 
+class TestMeshReadPropertiesTest(TestFixture.TestFixture):
+    def setUp(self):
+        self.props = ReadExodusMesh.read_exodus_mesh_element_properties(
+            pathlib.Path(__file__).parent.joinpath('read_material_property_test.exo'),
+            ['bulk', 'shear'], blockNum=1
+        )
+
+    def test_property_mins_and_maxs(self):
+        self.assertAlmostEqual(np.min(self.props, axis=0)[0], 0.26575326, 8)
+        self.assertAlmostEqual(np.min(self.props, axis=0)[1], 2.3917793, 8)
+        self.assertAlmostEqual(np.max(self.props, axis=0)[0], 1.38727616, 8)
+        self.assertAlmostEqual(np.max(self.props, axis=0)[1], 12.48548545, 8)
+
+    def test_bad_property_names(self):
+        with pytest.raises(KeyError):
+            self.props = ReadExodusMesh.read_exodus_mesh_element_properties(
+                pathlib.Path(__file__).parent.joinpath('read_material_property_test.exo'),
+                ['bulk1', 'shear1'], blockNum=1
+            )
+
 if __name__ == '__main__':
     unittest.main()