diff --git a/docs/source/io_formats/summary.rst b/docs/source/io_formats/summary.rst
index 1f30c69bb1..05c6fff9c0 100644
--- a/docs/source/io_formats/summary.rst
+++ b/docs/source/io_formats/summary.rst
@@ -34,7 +34,7 @@ The current version of the summary file format is 6.0.
              'material', 'universe', or 'lattice'.
            - **material** (*int* or *int[]*) -- Unique ID of the material(s)
              assigned to the cell. This dataset is present only if fill_type is
-             set to 'normal'.  The value '-1' signifies void material.  The data
+             set to 'normal'.  The value '-2' signifies void material.  The data
              is an array if the cell uses distributed materials, otherwise it is
              a scalar.
            - **temperature** (*double[]*) -- Temperature of the cell in Kelvin.
diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py
index 258e773dce..75986d2988 100644
--- a/openmc/capi/cell.py
+++ b/openmc/capi/cell.py
@@ -106,9 +106,10 @@ class Cell(_FortranObjectWithID):
 
         if fill_type.value == 1:
             if n.value > 1:
-                return [Material(index=i) for i in indices[:n.value]]
+                #TODO: off-by-one
+                return [Material(index=i+1) for i in indices[:n.value]]
             else:
-                return Material(index=indices[0])
+                return Material(index=indices[0]+1)
         else:
             raise NotImplementedError
 
@@ -116,14 +117,14 @@ class Cell(_FortranObjectWithID):
     def fill(self, fill):
         if isinstance(fill, Iterable):
             n = len(fill)
-            indices = (c_int32*n)(*(m._index if m is not None else -1
+            indices = (c_int32*n)(*(m._index if m is not None else -2
                                     for m in fill))
             _dll.openmc_cell_set_fill(self._index, 1, n, indices)
         elif isinstance(fill, Material):
             indices = (c_int32*1)(fill._index)
             _dll.openmc_cell_set_fill(self._index, 1, 1, indices)
         elif fill is None:
-            indices = (c_int32*1)(-1)
+            indices = (c_int32*1)(-2)
             _dll.openmc_cell_set_fill(self._index, 1, 1, indices)
 
     def set_temperature(self, T, instance=None):
diff --git a/src/api.F90 b/src/api.F90
index 6a6c10671a..60d94784de 100644
--- a/src/api.F90
+++ b/src/api.F90
@@ -37,8 +37,6 @@ module openmc_api
   public :: openmc_calculate_volumes
   public :: openmc_cell_filter_get_bins
   public :: openmc_cell_get_id
-  public :: openmc_cell_get_fill
-  public :: openmc_cell_set_fill
   public :: openmc_cell_set_id
   public :: openmc_cell_set_temperature
   public :: openmc_energy_filter_get_bins
diff --git a/src/cell.cpp b/src/cell.cpp
index 554ca74f1b..52b9a2e44f 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -8,6 +8,8 @@
 #include "error.h"
 #include "hdf5_interface.h"
 #include "lattice.h"
+#include "material.h"
+#include "openmc.h"
 #include "surface.h"
 #include "xml_interface.h"
 
@@ -219,10 +221,25 @@ Cell::Cell(pugi::xml_node cell_node)
     fill = C_NONE;
   }
 
+  // Read the material element.  There can be zero materials (filled with a
+  // universe), more than one material (distribmats), and some materials may
+  // be "void".
   if (check_for_node(cell_node, "material")) {
-    //TODO: read material ids.
-    material.push_back(C_NONE+1);
-    material.shrink_to_fit();
+    std::vector<std::string> mats
+         {get_node_array<std::string>(cell_node, "material")};
+    if (mats.size() > 0) {
+      material.reserve(mats.size());
+      for (std::string mat : mats) {
+        if (mat.compare("void") == 0) {
+          material.push_back(MATERIAL_VOID);
+        } else {
+          material.push_back(std::stoi(mat));
+        }
+      }
+    } else {
+      material.push_back(C_NONE);
+      material.shrink_to_fit();
+    }
   } else {
     material.push_back(C_NONE);
     material.shrink_to_fit();
@@ -460,6 +477,65 @@ read_cells(pugi::xml_node* node)
   }
 }
 
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+extern "C" int
+openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
+{
+  if (index >= 1 && index <= global_cells.size()) {
+    //TODO: off-by-one
+    Cell& c {*global_cells[index - 1]};
+    *type = c.type;
+    if (c.type == FILL_MATERIAL) {
+      *indices = c.material.data();
+      *n = c.material.size();
+    } else {
+      *indices = &c.fill;
+      *n = 1;
+    }
+  } else {
+    strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  return 0;
+}
+
+extern "C" int
+openmc_cell_set_fill(int32_t index, int type, int32_t n,
+                     const int32_t* indices)
+{
+  if (index >= 1 && index <= global_cells.size()) {
+    //TODO: off-by-one
+    Cell& c {*global_cells[index - 1]};
+    if (type == FILL_MATERIAL) {
+      c.type = FILL_MATERIAL;
+      c.material.clear();
+      for (int i = 0; i < n; i++) {
+        int i_mat = indices[i];
+        if (i_mat == MATERIAL_VOID) {
+          c.material.push_back(MATERIAL_VOID);
+        } else if (i_mat >= 1 && i_mat <= global_materials.size()) {
+          //TODO: off-by-one
+          c.material.push_back(i_mat - 1);
+        } else {
+          strcpy(openmc_err_msg, "Index in materials array is out of bounds.");
+          return OPENMC_E_OUT_OF_BOUNDS;
+        }
+      }
+    } else if (type == FILL_UNIVERSE) {
+      c.type = FILL_UNIVERSE;
+    } else {
+      c.type = FILL_LATTICE;
+    }
+  } else {
+    strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  return 0;
+}
+
 //==============================================================================
 // Fortran compatibility functions
 //==============================================================================
@@ -473,18 +549,23 @@ extern "C" {
 
   int cell_type(Cell* c) {return c->type;}
 
-  void cell_set_type(Cell* c, int type) {c->type = type;}
-
   int32_t cell_universe(Cell* c) {return c->universe;}
 
-  void cell_set_universe(Cell* c, int32_t universe) {c->universe = universe;}
-
   int32_t cell_fill(Cell* c) {return c->fill;}
 
-  int32_t* cell_fill_ptr(Cell* c) {return &c->fill;}
-
   int32_t cell_n_instances(Cell* c) {return c->n_instances;}
 
+  int cell_material_size(Cell* c) {return c->material.size();}
+
+  //TODO: off-by-one
+  int32_t cell_material(Cell* c, int i)
+  {
+    int32_t mat = c->material[i-1];
+    if (mat == C_NONE) return F90_NONE;
+    if (mat == MATERIAL_VOID) return MATERIAL_VOID;
+    return mat + 1;
+  }
+
   bool cell_simple(Cell* c) {return c->simple;}
 
   bool cell_contains(Cell* c, double xyz[3], double uvw[3], int32_t on_surface)
diff --git a/src/constants.F90 b/src/constants.F90
index ab85b277a6..4b0f392fd7 100644
--- a/src/constants.F90
+++ b/src/constants.F90
@@ -129,7 +129,7 @@ module constants
   integer(C_INT), bind(C, name='FILL_LATTICE')  :: FILL_LATTICE_C  = FILL_LATTICE
 
   ! Void material
-  integer, parameter :: MATERIAL_VOID = -1
+  integer, parameter :: MATERIAL_VOID = -2
 
   ! Flag to say that the outside of a lattice is not defined
   integer, parameter :: NO_OUTER_UNIVERSE = -1
diff --git a/src/constants.h b/src/constants.h
index b6bb82d375..72c72bad3f 100644
--- a/src/constants.h
+++ b/src/constants.h
@@ -112,7 +112,7 @@ constexpr char SUBSHELLS[][4]  {
 
 // Void material
 // TODO: refactor and remove
-constexpr int MATERIAL_VOID {-1};
+constexpr int MATERIAL_VOID {-2};
 
 // ============================================================================
 // CROSS SECTION RELATED CONSTANTS
@@ -434,6 +434,7 @@ constexpr int DIFF_NUCLIDE_DENSITY {2};
 constexpr int DIFF_TEMPERATURE {3};
 
 constexpr int C_NONE {-1};
+constexpr int F90_NONE {0};
 
 // Interpolation rules
 enum class Interpolation {
diff --git a/src/geometry.F90 b/src/geometry.F90
index 5bf0012a38..81e123d174 100644
--- a/src/geometry.F90
+++ b/src/geometry.F90
@@ -172,7 +172,7 @@ contains
           p % last_sqrtkT = p % sqrtkT
 
           ! Get distributed offset
-          if (size(c % material) > 1 .or. size(c % sqrtkT) > 1) then
+          if (c % material_size() > 1 .or. size(c % sqrtkT) > 1) then
             ! Distributed instances of this cell have different
             ! materials/temperatures. Determine which instance this is for
             ! assigning the matching material/temperature.
@@ -204,7 +204,7 @@ contains
           end if
 
           ! Save the material
-          if (size(c % material) > 1) then
+          if (c % material_size() > 1) then
             p % material = c % material(offset + 1)
           else
             p % material = c % material(1)
diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp
index 3f36cf299a..8399817fa4 100644
--- a/src/geometry_aux.cpp
+++ b/src/geometry_aux.cpp
@@ -8,6 +8,7 @@
 #include "constants.h"
 #include "error.h"
 #include "lattice.h"
+#include "material.h"
 
 
 namespace openmc {
@@ -15,7 +16,7 @@ namespace openmc {
 //==============================================================================
 
 void
-adjust_indices_c()
+adjust_indices()
 {
   // Adjust material/fill idices.
   for (Cell *c : global_cells) {
@@ -36,8 +37,21 @@ adjust_indices_c()
         fatal_error(err_msg);
       }
     } else {
-      //TODO: materials
       c->type = FILL_MATERIAL;
+      for (auto it = c->material.begin(); it != c->material.end(); it++) {
+        int32_t mid = *it;
+        if (mid != MATERIAL_VOID) {
+          auto search = material_map.find(mid);
+          if (search != material_map.end()) {
+            *it = search->second;
+          } else {
+            std::stringstream err_msg;
+            err_msg << "Could not find material " << mid
+                    << " specified on cell " << c->id;
+            fatal_error(err_msg);
+          }
+        }
+      }
     }
   }
 
diff --git a/src/geometry_aux.h b/src/geometry_aux.h
index b998ad32c9..fb5940d1c9 100644
--- a/src/geometry_aux.h
+++ b/src/geometry_aux.h
@@ -13,7 +13,7 @@ namespace openmc {
 //! Replace Universe, Lattice, and Material IDs with indices.
 //==============================================================================
 
-extern "C" void adjust_indices_c();
+extern "C" void adjust_indices();
 
 //==============================================================================
 //! Figure out which Universe is the root universe.
diff --git a/src/geometry_header.F90 b/src/geometry_header.F90
index 5a2dca282f..93bc9918b2 100644
--- a/src/geometry_header.F90
+++ b/src/geometry_header.F90
@@ -40,12 +40,6 @@ module geometry_header
       integer(C_INT)                 :: type
     end function cell_type_c
 
-    subroutine cell_set_type_c(cell_ptr, type) bind(C, name='cell_set_type')
-      import C_PTR, C_INT
-      type(C_PTR),    intent(in), value :: cell_ptr
-      integer(C_INT), intent(in), value :: type
-    end subroutine cell_set_type_c
-
     function cell_universe_c(cell_ptr) bind(C, name='cell_universe') &
          result(universe)
       import C_PTR, C_INT32_T
@@ -53,25 +47,12 @@ module geometry_header
       integer(C_INT32_T)             :: universe
     end function cell_universe_c
 
-    subroutine cell_set_universe_c(cell_ptr, universe) &
-         bind(C, name='cell_set_universe')
-      import C_PTR, C_INT32_T
-      type(C_PTR),        intent(in), value :: cell_ptr
-      integer(C_INT32_T), intent(in), value :: universe
-    end subroutine cell_set_universe_c
-
     function cell_fill_c(cell_ptr) bind(C, name="cell_fill") result(fill)
       import C_PTR, C_INT32_T
       type(C_PTR), intent(in), value :: cell_ptr
       integer(C_INT32_T)             :: fill
     end function cell_fill_c
 
-    function cell_fill_ptr(cell_ptr) bind(C) result(fill_ptr)
-      import C_PTR
-      type(C_PTR), intent(in), value :: cell_ptr
-      type(C_PTR)                    :: fill_ptr
-    end function cell_fill_ptr
-
     function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') &
          result(n_instances)
       import C_PTR, C_INT32_T
@@ -79,6 +60,21 @@ module geometry_header
       integer(C_INT32_T)             :: n_instances
     end function cell_n_instances_c
 
+    function cell_material_size_c(cell_ptr) bind(C, name='cell_material_size') &
+         result(n)
+      import C_PTR, C_INT
+      type(C_PTR), intent(in), value :: cell_ptr
+      integer(C_INT)                 :: n
+    end function cell_material_size_c
+
+    function cell_material_c(cell_ptr, i) bind(C, name='cell_material') &
+         result(mat)
+      import C_PTR, C_INT, C_INT32_T
+      type(C_PTR),    intent(in), value :: cell_ptr
+      integer(C_INT), intent(in), value :: i
+      integer(C_INT32_T)                :: mat
+    end function cell_material_c
+
     function cell_simple_c(cell_ptr) bind(C, name='cell_simple') result(simple)
       import C_PTR, C_BOOL
       type(C_PTR), intent(in), value :: cell_ptr
@@ -253,9 +249,6 @@ module geometry_header
   type Cell
     type(C_PTR) :: ptr
 
-    integer, allocatable :: material(:)    ! Material within cell.  Multiple
-                                           !  materials for distribcell
-                                           !  instances.  0 signifies a universe
     integer, allocatable :: region(:)      ! Definition of spatial region as
                                            !  Boolean expression of half-spaces
     integer :: distribcell_index           ! Index corresponding to this cell in
@@ -274,11 +267,11 @@ module geometry_header
     procedure :: id => cell_id
     procedure :: set_id => cell_set_id
     procedure :: type => cell_type
-    procedure :: set_type => cell_set_type
     procedure :: universe => cell_universe
-    procedure :: set_universe => cell_set_universe
     procedure :: fill => cell_fill
     procedure :: n_instances => cell_n_instances
+    procedure :: material_size => cell_material_size
+    procedure :: material => cell_material
     procedure :: simple => cell_simple
     procedure :: distance => cell_distance
     procedure :: offset => cell_offset
@@ -389,24 +382,12 @@ contains
     type = cell_type_c(this % ptr)
   end function cell_type
 
-  subroutine cell_set_type(this, type)
-    class(Cell),    intent(in) :: this
-    integer(C_INT), intent(in) :: type
-    call cell_set_type_c(this % ptr, type)
-  end subroutine cell_set_type
-
   function cell_universe(this) result(universe)
     class(Cell), intent(in) :: this
     integer(C_INT32_T)      :: universe
     universe = cell_universe_c(this % ptr)
   end function cell_universe
 
-  subroutine cell_set_universe(this, universe)
-    class(Cell),        intent(in) :: this
-    integer(C_INT32_T), intent(in) :: universe
-    call cell_set_universe_c(this % ptr, universe)
-  end subroutine cell_set_universe
-
   function cell_fill(this) result(fill)
     class(Cell), intent(in) :: this
     integer(C_INT32_T)      :: fill
@@ -419,6 +400,19 @@ contains
     n_instances = cell_n_instances_c(this % ptr)
   end function cell_n_instances
 
+  function cell_material_size(this) result(n)
+    class(Cell), intent(in) :: this
+    integer(C_INT)          :: n
+    n = cell_material_size_c(this % ptr)
+  end function cell_material_size
+
+  function cell_material(this, i) result(mat)
+    class(Cell), intent(in) :: this
+    integer,     intent(in) :: i
+    integer(C_INT32_T)      :: mat
+    mat = cell_material_c(this % ptr, i)
+  end function cell_material
+
   function cell_simple(this) result(simple)
     class(Cell), intent(in) :: this
     logical(C_BOOL)         :: simple
@@ -468,7 +462,7 @@ contains
     if (present(sab_temps)) allocate(sab_temps(n_sab_tables))
 
     do i = 1, size(cells)
-      do j = 1, size(cells(i) % material)
+      do j = 1, cells(i) % material_size()
         ! Skip any non-material cells and void materials
         if (cells(i) % material(j) == NONE .or. &
              cells(i) % material(j) == MATERIAL_VOID) cycle
@@ -598,33 +592,6 @@ contains
   end function openmc_get_cell_index
 
 
-  function openmc_cell_get_fill(index, type, indices, n) result(err) bind(C)
-    integer(C_INT32_T), value, intent(in) :: index
-    integer(C_INT), intent(out) :: type
-    integer(C_INT32_T), intent(out) :: n
-    type(C_PTR), intent(out) :: indices
-    integer(C_INT) :: err
-
-    err = 0
-    if (index >= 1 .and. index <= size(cells)) then
-      associate (c => cells(index))
-        type = c % type()
-        select case (type)
-        case (FILL_MATERIAL)
-          n = size(c % material)
-          indices = C_LOC(c % material(1))
-        case (FILL_UNIVERSE, FILL_LATTICE)
-          n = 1
-          indices = cell_fill_ptr(c % ptr)
-        end select
-      end associate
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg("Index in cells array is out of bounds.")
-    end if
-  end function openmc_cell_get_fill
-
-
   function openmc_cell_get_id(index, id) result(err) bind(C)
     ! Return the ID of a cell
     integer(C_INT32_T), value       :: index
@@ -641,49 +608,6 @@ contains
   end function openmc_cell_get_id
 
 
-  function openmc_cell_set_fill(index, type, n, indices) result(err) bind(C)
-    ! Set the fill for a cell
-    integer(C_INT32_T), value, intent(in) :: index    ! index in cells
-    integer(C_INT), value, intent(in)     :: type
-    integer(c_INT32_T), value, intent(in) :: n
-    integer(C_INT32_T), intent(in) :: indices(n)
-    integer(C_INT) :: err
-
-    integer :: i, j
-
-    err = 0
-    if (index >= 1 .and. index <= size(cells)) then
-      associate (c => cells(index))
-        select case (type)
-        case (FILL_MATERIAL)
-          if (allocated(c % material)) deallocate(c % material)
-          allocate(c % material(n))
-
-          call c % set_type(FILL_MATERIAL)
-          do i = 1, n
-            j = indices(i)
-            if ((j >= 1 .and. j <= n_materials) .or. j == MATERIAL_VOID) then
-              c % material(i) = j
-            else
-              err = E_OUT_OF_BOUNDS
-              call set_errmsg("Index " // trim(to_str(j)) // " in the &
-                   &materials array is out of bounds.")
-            end if
-          end do
-        case (FILL_UNIVERSE)
-          call c % set_type(FILL_UNIVERSE)
-        case (FILL_LATTICE)
-          call c % set_type(FILL_LATTICE)
-        end select
-      end associate
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg("Index in cells array is out of bounds.")
-    end if
-
-  end function openmc_cell_set_fill
-
-
   function openmc_cell_set_id(index, id) result(err) bind(C)
     ! Set the ID of a cell
     integer(C_INT32_T), value, intent(in) :: index
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index aef9a14c07..c10d31deab 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -48,8 +48,8 @@ module input_xml
   save
 
   interface
-    subroutine adjust_indices_c() bind(C)
-    end subroutine adjust_indices_c
+    subroutine adjust_indices() bind(C)
+    end subroutine adjust_indices
 
     subroutine allocate_offset_tables(n_maps) bind(C)
       import C_INT
@@ -1114,42 +1114,6 @@ contains
              // to_str(c % id()))
       end if
 
-      ! Read material
-      if (check_for_node(node_cell, "material")) then
-        n_mats = node_word_count(node_cell, "material")
-
-        if (n_mats > 0) then
-          allocate(sarray(n_mats))
-          call get_node_array(node_cell, "material", sarray)
-
-          allocate(c % material(n_mats))
-          do j = 1, n_mats
-            select case(trim(to_lower(sarray(j))))
-            case ('void')
-              c % material(j) = MATERIAL_VOID
-            case default
-              c % material(j) = int(str_to_int(sarray(j)), 4)
-
-              ! Check for error
-              if (c % material(j) == ERROR_INT) then
-                call fatal_error("Invalid material specified on cell " &
-                     // to_str(c % id()))
-              end if
-            end select
-          end do
-
-          deallocate(sarray)
-
-        else
-          allocate(c % material(1))
-          c % material(1) = NONE
-        end if
-
-      else
-        allocate(c % material(1))
-        c % material(1) = NONE
-      end if
-
       ! Check for region specification (also under deprecated name surfaces)
       if (check_for_node(node_cell, "surfaces")) then
         call warning("The use of 'surfaces' is deprecated and will be &
@@ -1550,7 +1514,7 @@ contains
     ! Get pointer to list of XML <material>
     call get_node_list(root, "material", node_mat_list)
 
-    ! Allocate cells array
+    ! Allocate materials array
     n_materials = size(node_mat_list)
     allocate(materials(n_materials))
     allocate(material_temps(n_materials))
@@ -1575,12 +1539,6 @@ contains
              mat % depletable = .true.
       end if
 
-      ! Check to make sure 'id' hasn't been used
-      if (material_dict % has(mat % id())) then
-        call fatal_error("Two or more materials use the same unique ID: " &
-             // to_str(mat % id()))
-      end if
-
       ! Copy material name
       if (check_for_node(node_mat, "name")) then
         call get_node_value(node_mat, "name", mat % name)
@@ -3861,10 +3819,10 @@ contains
 
       ! Set the number of temperatures equal to the number of materials.
       deallocate(cells(i) % sqrtkT)
-      allocate(cells(i) % sqrtkT(size(cells(i) % material)))
+      allocate(cells(i) % sqrtkT(cells(i) % material_size()))
 
       ! Check each of the cell materials for temperature data.
-      do j = 1, size(cells(i) % material)
+      do j = 1, cells(i) % material_size()
         ! Arbitrarily set void regions to 0K.
         if (cells(i) % material(j) == MATERIAL_VOID) then
           cells(i) % sqrtkT(j) = ZERO
@@ -3932,45 +3890,6 @@ contains
 
   end subroutine read_multipole_data
 
-!===============================================================================
-! ADJUST_INDICES changes the values for 'surfaces' for each cell and the
-! material index assigned to each to the indices in the surfaces and material
-! array rather than the unique IDs assigned to each surface and material. Also
-! assigns boundary conditions to surfaces based on those read into the bc_dict
-! dictionary
-!===============================================================================
-
-  subroutine adjust_indices()
-
-    integer :: i                      ! index for various purposes
-    integer :: j                      ! index for various purposes
-    integer :: id                     ! user-specified id
-
-    call adjust_indices_c()
-
-    do i = 1, n_cells
-      associate (c => cells(i))
-
-      ! =======================================================================
-      ! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL
-
-      if (c % material(1) /= NONE) then
-        do j = 1, size(c % material)
-          id = c % material(j)
-          if (id == MATERIAL_VOID) then
-          else if (material_dict % has(id)) then
-            c % material(j) = material_dict % get(id)
-          else
-            call fatal_error("Could not find material " // trim(to_str(id)) &
-                 // " specified on cell " // trim(to_str(c % id())))
-          end if
-        end do
-      end if
-      end associate
-    end do
-
-  end subroutine adjust_indices
-
 !===============================================================================
 ! PREPARE_DISTRIBCELL initializes any distribcell filters present and sets the
 ! offsets for distribcells
@@ -3994,7 +3913,7 @@ contains
 
     ! Find all cells with multiple (distributed) materials or temperatures.
     do i = 1, n_cells
-      if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then
+      if (cells(i) % material_size() > 1 .or. size(cells(i) % sqrtkT) > 1) then
         call cell_list % add(i)
       end if
     end do
@@ -4003,10 +3922,10 @@ contains
     ! number of respective cell instances.
     do i = 1, n_cells
       associate (c => cells(i))
-        if (size(c % material) > 1) then
-          if (size(c % material) /= c % n_instances()) then
+        if (c % material_size() > 1) then
+          if (c % material_size() /= c % n_instances()) then
             call fatal_error("Cell " // trim(to_str(c % id())) // " was &
-                 &specified with " // trim(to_str(size(c % material))) &
+                 &specified with " // trim(to_str(c % material_size())) &
                  // " materials but has " // trim(to_str(c % n_instances())) &
                  // " distributed instances. The number of materials must &
                  &equal one or the number of instances.")
diff --git a/src/material.cpp b/src/material.cpp
index 49f0e24158..3bc2452d39 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -1,5 +1,8 @@
 #include "material.h"
 
+#include <string>
+#include <sstream>
+
 #include "error.h"
 #include "xml_interface.h"
 
@@ -14,6 +17,7 @@ namespace openmc {
 //==============================================================================
 
 std::vector<Material*> global_materials;
+std::unordered_map<int32_t, int32_t> material_map;
 
 //==============================================================================
 // Material implementation
@@ -26,8 +30,6 @@ Material::Material(pugi::xml_node material_node)
   } else {
     fatal_error("Must specify id of material in materials XML file.");
   }
-
-  std::cout << "Reading material with id " << id << std::endl;
 }
 
 //==============================================================================
@@ -41,6 +43,19 @@ read_materials(pugi::xml_node* node)
   for (pugi::xml_node material_node: node->children("material")) {
     global_materials.push_back(new Material(material_node));
   }
+
+  // Populate the material map.
+  for (int i = 0; i < global_materials.size(); i++) {
+    int32_t mid = global_materials[i]->id;
+    auto search = material_map.find(mid);
+    if (search == material_map.end()) {
+      material_map[mid] = i;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Two or more materials use the same unique ID: " << mid;
+      fatal_error(err_msg);
+    }
+  }
 }
 
 //==============================================================================
@@ -61,6 +76,13 @@ extern "C" {
       global_materials.push_back(new Material());
     }
   }
+
+  void free_memory_material_c()
+  {
+    for (Material *mat : global_materials) {delete mat;}
+    global_materials.clear();
+    material_map.clear();
+  }
 }
 
 } // namespace openmc
diff --git a/src/material.h b/src/material.h
index f2e3df2fb4..0a9a099a44 100644
--- a/src/material.h
+++ b/src/material.h
@@ -1,6 +1,7 @@
-#ifndef OPENMC_CELL_H
-#define OPENMC_CELL_H
+#ifndef OPENMC_MATERIAL_H
+#define OPENMC_MATERIAL_H
 
+#include <unordered_map>
 #include <vector>
 
 #include "pugixml.hpp"
@@ -14,6 +15,7 @@ namespace openmc {
 
 class Material;
 extern std::vector<Material*> global_materials;
+extern std::unordered_map<int32_t, int32_t> material_map;
 
 //==============================================================================
 //! A substance with constituent nuclides and thermal scattering data
@@ -30,4 +32,4 @@ public:
 };
 
 } // namespace openmc
-#endif // OPENMC_CELL_H
+#endif // OPENMC_MATERIAL_H
diff --git a/src/material_header.F90 b/src/material_header.F90
index 5bfbadfcf9..3a68a48694 100644
--- a/src/material_header.F90
+++ b/src/material_header.F90
@@ -484,6 +484,11 @@ contains
 !===============================================================================
 
   subroutine free_memory_material()
+    interface
+      subroutine free_memory_material_c() bind(C)
+      end subroutine free_memory_material_c
+    end interface
+    call free_memory_material_c()
     n_materials = 0
     if (allocated(materials)) deallocate(materials)
     call material_dict % clear()
diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90
index 1b3ed35113..f6289518d2 100644
--- a/src/mgxs_data.F90
+++ b/src/mgxs_data.F90
@@ -151,7 +151,7 @@ contains
     allocate(kTs(size(materials)))
 
     do i = 1, size(cells)
-      do j = 1, size(cells(i) % material)
+      do j = 1, cells(i) % material_size()
 
         ! Skip any non-material cells and void materials
         if (cells(i) % material(j) == NONE .or. &
diff --git a/src/summary.F90 b/src/summary.F90
index 507642f713..13bc426637 100644
--- a/src/summary.F90
+++ b/src/summary.F90
@@ -198,7 +198,7 @@ contains
       case (FILL_MATERIAL)
         call write_dataset(cell_group, "fill_type", "material")
 
-        if (size(c % material) == 1) then
+        if (c % material_size() == 1) then
           if (c % material(1) == MATERIAL_VOID) then
             call write_dataset(cell_group, "material", MATERIAL_VOID)
           else
@@ -206,8 +206,8 @@ contains
                  materials(c % material(1)) % id())
           end if
         else
-          allocate(cell_materials(size(c % material)))
-          do j = 1, size(c % material)
+          allocate(cell_materials(c % material_size()))
+          do j = 1, c % material_size()
             if (c % material(j) == MATERIAL_VOID) then
               cell_materials(j) = MATERIAL_VOID
             else
diff --git a/src/xml_interface.h b/src/xml_interface.h
index f10143697a..823948de2f 100644
--- a/src/xml_interface.h
+++ b/src/xml_interface.h
@@ -17,7 +17,7 @@ check_for_node(pugi::xml_node node, const char *name)
 }
 
 std::string get_node_value(pugi::xml_node node, const char *name,
-  bool lowercase=false, bool strip=false);
+                           bool lowercase=true, bool strip=true);
 
 template <typename T>
 std::vector<T> get_node_array(pugi::xml_node node, const char* name)