diff --git a/n88tools/pistoia.py b/n88tools/pistoia.py
index 35e4f16..3806210 100644
--- a/n88tools/pistoia.py
+++ b/n88tools/pistoia.py
@@ -10,7 +10,7 @@
 
 from __future__ import division
 import sys
-from .N88ReportedError import N88ReportedError
+from n88tools.N88ReportedError import N88ReportedError
 import numpy
 from numpy.core import *
 
@@ -249,7 +249,7 @@ def StatsTable (data, stats, labels=None):
             values = []
             for i in range(num_col):
                 values.append (CalculateStats(data[:,i], stats))
-            for v in values[0].keys():
+            for v in list(values[0].keys()):
                 row_values = []
                 for i in range(num_col):
                     row_values.append (values[i][v])
@@ -335,28 +335,37 @@ def StatsTable (data, stats, labels=None):
     matNames = materialTable.variables['MaterialName'][:]
     for id,matName in zip (materialIds, matNames):
         material = materialDefinitions.groups[matName]
-        if material.Type == "LinearIsotropic":
-            modulus_value = material.E
-        elif material.Type == "LinearOrthotropic":
-            modulus_value = max(material.E)
+        if hasattr(material, 'E'):
+            if material.Type == "LinearIsotropic":
+                modulus_value = material.E
+            elif material.Type == "LinearOrthotropic":
+                modulus_value = max(material.E)
+            else:
+                raise N88ReportedError ("ERROR: Unsupported material type for Pistoia calculation: %s" % material.Type)
+            indices = numpy.nonzero(elementMaterials == id)
+            modulus[indices] = modulus_value
         else:
-            raise N88ReportedError ("ERROR: Unsupported material type for Pistoia calculation: %s" % material.Type)
-        indices = numpy.nonzero(elementMaterials == id)
-        modulus[indices] = modulus_value
+            if material.Type != "LinearIsotropic":
+                raise N88ReportedError ("ERROR: Unsupported material type for Pistoia calculation: %s" % material.Type)
+
+            modulus_values = numpy.array(material.variables['E'][:])
+            for i, E in enumerate(modulus_values):
+                indices = numpy.nonzero(elementMaterials == i)
+                modulus[indices] = E            
     sed = zeros (elementValues.variables["StrainEnergyDensity"].shape, float64)
     sed[:] = elementValues.variables["StrainEnergyDensity"][:]
-    if mask != None:
+    if mask is not None:
         sed = numpy.ma.MaskedArray (sed, mask)
         modulus = numpy.ma.MaskedArray (modulus, mask)
     ees = sqrt(2*sed/modulus)
     del sed
-    if mask == None:
+    if mask is None:
         numberOfSelectedElements = numberOfElements
     else:
         numberOfSelectedElements = ees.count()
     if numberOfSelectedElements == 0:
         raise N88ReportedError ("ERROR: No elements in selection.")
-    if mask == None:
+    if mask is None:
         sort_indices = argsort (ees)
     else:
         sort_indices = ees.argsort(fill_value=inf)[:numberOfSelectedElements]
@@ -419,9 +428,11 @@ def StatsTable (data, stats, labels=None):
     out.write ("Distribution of failed materials.\n")
     out.write (subtable_delimiter)
     g_mat_failed = elementMaterials[sort_indices[el:]]
-    nels_failed = zeros (sizeOfMaterialTable, int)
-    for m in range(sizeOfMaterialTable):
-        nels_failed[m] = sum (g_mat_failed == materialIds[m])
+    unique_materials = numpy.unique(elementMaterials)
+    nMaterials = len(unique_materials)
+    nels_failed = zeros (nMaterials, int)
+    for i,m in enumerate(unique_materials):
+        nels_failed[i] = sum(g_mat_failed == m)
     total = sum(nels_failed)
     col_width = 12
     left_margin = (width - 3*col_width)//2
@@ -429,8 +440,8 @@ def StatsTable (data, stats, labels=None):
     table_line = "-"*10
     row_format = " "*left_margin + "%%%ds" % col_width + "%%%dd" % col_width + "%%%d.2f" % col_width + "\n"
     out.write (heading_format % ("material","# els","%"))
-    for m in range(sizeOfMaterialTable):
-        out.write (row_format % (materialIds[m], nels_failed[m], numpy.float64(100)*nels_failed[m]/total))
+    for i,m in enumerate(unique_materials):
+        out.write (row_format % (m, nels_failed[i], numpy.float64(100)*nels_failed[i]/total))
     out.write (heading_format % ((table_line,)*3))
     out.write (row_format % ("total",total,100))