CarpetLib: Use AT to index std::vector

Define a macro AT() to index into std::vector.  Depending on the macro
NDEBUG, AT() is defined either as at(), providing index checking, or
as operator[], providing no checking.

darcs-hash:20070203205854-dae7b-a1999c88c95ba12b1ee66505f712aefdd67d7e6f.gz

14 files changed, 220 insertions, 208 deletions

diff --git a/Carpet/CarpetLib/src/commstate.cc b/Carpet/CarpetLib/src/commstate.cc
index 0ac7b21de..9eecec8f2 100644
--- a/Carpet/CarpetLib/src/commstate.cc
+++ b/Carpet/CarpetLib/src/commstate.cc
@@ -48,8 +48,8 @@ comm_state::comm_state ()
   {                                                                    \
     T dummy;                                                           \
     int type = dist::c_datatype (dummy);                               \
-    typebufs.at(type).datatypesize = sizeof (dummy);                   \
-    typebufs.at(type).mpi_datatype = dist::datatype (dummy);           \
+    typebufs.AT(type).datatypesize = sizeof (dummy);                   \
+    typebufs.AT(type).mpi_datatype = dist::datatype (dummy);           \
   }
 #include "instantiate"
 #undef INSTANTIATE
@@ -253,15 +253,15 @@ bool comm_state::AllPostedCommunicationsFinished ()
     if (reduce_mpi_waitall) {
       size_t nreqs = 0;
       for (size_t i=0; i<srequests.size(); ++i) {
-        if (srequests.at(i) != MPI_REQUEST_NULL) {
+        if (srequests.AT(i) != MPI_REQUEST_NULL) {
           ++nreqs;
         }
       }
       vector<MPI_Request> reqs(nreqs);
       nreqs = 0;
       for (size_t i=0; i<srequests.size(); ++i) {
-        if (srequests.at(i) != MPI_REQUEST_NULL) {
-          reqs.at(nreqs) = srequests.at(i);
+        if (srequests.AT(i) != MPI_REQUEST_NULL) {
+          reqs.AT(nreqs) = srequests.AT(i);
           ++nreqs;
         }
       }
@@ -284,22 +284,22 @@ bool comm_state::AllPostedCommunicationsFinished ()
   if (use_waitall) {
     // mark all posted recveive buffers as ready
     for (size_t i = 0; i < recvbuffers_ready.size(); i++) {
-      recvbuffers_ready.at(i) = rrequests.at(i) != MPI_REQUEST_NULL;
+      recvbuffers_ready.AT(i) = rrequests.AT(i) != MPI_REQUEST_NULL;
     }
 
     // wait for completion of all posted receive operations
     if (reduce_mpi_waitall) {
       size_t nreqs = 0;
       for (size_t i=0; i<rrequests.size(); ++i) {
-        if (rrequests.at(i) != MPI_REQUEST_NULL) {
+        if (rrequests.AT(i) != MPI_REQUEST_NULL) {
           ++nreqs;
         }
       }
       vector<MPI_Request> reqs(nreqs);
       nreqs = 0;
       for (size_t i=0; i<rrequests.size(); ++i) {
-        if (rrequests.at(i) != MPI_REQUEST_NULL) {
-          reqs.at(nreqs) = rrequests.at(i);
+        if (rrequests.AT(i) != MPI_REQUEST_NULL) {
+          reqs.AT(nreqs) = rrequests.AT(i);
           ++nreqs;
         }
       }
@@ -327,8 +327,8 @@ bool comm_state::AllPostedCommunicationsFinished ()
 
     // mark the recveive buffers of completed communications as ready
     for (int i = 0; i < num_completed_recvs_; i++) {
-      assert (rrequests.at(completed_recvs.at(i)) == MPI_REQUEST_NULL);
-      recvbuffers_ready.at(completed_recvs.at(i)) = true;
+      assert (rrequests.AT(completed_recvs.AT(i)) == MPI_REQUEST_NULL);
+      recvbuffers_ready.AT(completed_recvs.AT(i)) = true;
     }
   }
 
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 9067d2942..e3feae50f 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -453,13 +453,13 @@ void data <T>
   
   vector <data const *> srcs (gsrcs.size());
   for (size_t t=0; t<srcs.size(); ++t) {
-    srcs.at(t) = dynamic_cast <data const *> (gsrcs.at(t));
+    srcs.AT(t) = dynamic_cast <data const *> (gsrcs.AT(t));
   }
   assert (srcs.size() == times.size() and srcs.size() > 0);
   
   for (size_t t=0; t<srcs.size(); ++t) {
-    assert (srcs.at(t)->has_storage());
-    assert (proc() == srcs.at(t)->proc());
+    assert (srcs.AT(t)->has_storage());
+    assert (proc() == srcs.AT(t)->proc());
   }
   
   assert (dist::rank() == proc());
@@ -508,7 +508,7 @@ void data <T>
   }
   if (timelevel == -1) {
     for (size_t tl=0; tl<times.size(); ++tl) {
-      if (abs (times.at(tl) - time) < eps) {
+      if (abs (times.AT(tl) - time) < eps) {
         timelevel = tl;
         break;
       }
@@ -522,24 +522,24 @@ void data <T>
     
     for (size_t tl=0; tl<times.size(); ++tl) {
       
-      tmps.at(tl) =
+      tmps.AT(tl) =
         new data (this->varindex, this->cent, this->transport_operator);
-      tmps.at(tl)->allocate (box, this->proc());
+      tmps.AT(tl)->allocate (box, this->proc());
       
-      tmps.at(tl)->interpolate_p_r (srcs.at(tl), box, order_space);
+      tmps.AT(tl)->interpolate_p_r (srcs.AT(tl), box, order_space);
       
     }
     
     time_interpolate (tmps, box, times, time, order_time);
     
     for (size_t tl=0; tl<times.size(); ++tl) {
-      delete tmps.at(tl);
+      delete tmps.AT(tl);
     }
     
   } else {
     // No time interpolation
     
-    interpolate_p_r (srcs.at(timelevel), box, order_space);
+    interpolate_p_r (srcs.AT(timelevel), box, order_space);
     
   } // if
 }
@@ -848,32 +848,32 @@ void data <T>
     
   case 1:
     assert (times.size() >= 2);
-    interpolate_3d_2tl (static_cast <T const *> (srcs.at(0)->storage()),
-                        times.at(0),
-                        static_cast <T const *> (srcs.at(1)->storage()),
-                        times.at(1),
-                        srcs.at(0)->shape(),
+    interpolate_3d_2tl (static_cast <T const *> (srcs.AT(0)->storage()),
+                        times.AT(0),
+                        static_cast <T const *> (srcs.AT(1)->storage()),
+                        times.AT(1),
+                        srcs.AT(0)->shape(),
                         static_cast <T *> (this->storage()),
                         time,
                         this->shape(),
-                        srcs.at(0)->extent(),
+                        srcs.AT(0)->extent(),
                         this->extent(),
                         box);
     break;
     
   case 2:
     assert (times.size() >= 3);
-    interpolate_3d_3tl (static_cast <T const *> (srcs.at(0)->storage()),
-                        times.at(0),
-                        static_cast <T const *> (srcs.at(1)->storage()),
-                        times.at(1),
-                        static_cast <T const *> (srcs.at(2)->storage()),
-                        times.at(2),
-                        srcs.at(0)->shape(),
+    interpolate_3d_3tl (static_cast <T const *> (srcs.AT(0)->storage()),
+                        times.AT(0),
+                        static_cast <T const *> (srcs.AT(1)->storage()),
+                        times.AT(1),
+                        static_cast <T const *> (srcs.AT(2)->storage()),
+                        times.AT(2),
+                        srcs.AT(0)->shape(),
                         static_cast <T *> (this->storage()),
                         time,
                         this->shape(),
-                        srcs.at(0)->extent(),
+                        srcs.AT(0)->extent(),
                         this->extent(),
                         box);
     break;
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index 3ae0b8261..afbabfba5 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -36,6 +36,17 @@ using namespace std;
 
 
 
+// Use this macro AT instead of vector's operator[] or at().
+// Depending on the macro NDEBUG, this macro AT either checks for
+// valid indices or not.
+#ifndef NDEBUG
+#  define AT(index) at(index)
+#else
+#  define AT(index) operator[](index)
+#endif
+
+
+
 // Begin a new line without flushing the output buffer
 char const * const eol = "\n";
 
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 67c2a94f7..6fdde0a20 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -108,12 +108,12 @@ void dh::allocate_bboxes ()
 {
   boxes.resize(h.mglevels());
   for (int ml=0; ml<h.mglevels(); ++ml) {
-    boxes.at(ml).resize(h.reflevels());
+    boxes.AT(ml).resize(h.reflevels());
     for (int rl=0; rl<h.reflevels(); ++rl) {
-      boxes.at(ml).at(rl).resize(h.components(rl));
+      boxes.AT(ml).AT(rl).resize(h.components(rl));
       for (int c=0; c<h.components(rl); ++c) {
         const ibbox intr = h.extent(ml,rl,c);
-        dboxes & b = boxes.at(ml).at(rl).at(c);
+        dboxes & b = boxes.AT(ml).AT(rl).AT(c);
 
 
 
@@ -132,7 +132,7 @@ void dh::allocate_bboxes ()
           for (int d=0; d<dim; ++d) {
             if (h.outer_boundaries(rl,c)[f][d]) {
               dist[f][d] = 0;
-              boxes.at(ml).at(rl).at(c).is_interproc[d][f] = false;
+              boxes.AT(ml).AT(rl).AT(c).is_interproc[d][f] = false;
             } else {
               bool const is_empty = intr.lower()[d] > intr.upper()[d];
               if (is_empty) {
@@ -147,7 +147,7 @@ void dh::allocate_bboxes ()
                 bnd -= h.extent(ml,rl,cc);
               }
               bool const is_interproc = bnd.empty();
-              boxes.at(ml).at(rl).at(c).is_interproc[d][f] = is_interproc;
+              boxes.AT(ml).AT(rl).AT(c).is_interproc[d][f] = is_interproc;
               if (! is_empty and ! is_interproc) {
                 dist[f][d] += buffers[f][d];
                 odist[f][d] = dist[f][d];
@@ -171,19 +171,19 @@ void dh::allocate_bboxes ()
       // Make owned regions disjoint
       for (int c=0; c<h.components(rl); ++c) {
         const ibbox intr = h.extent(ml,rl,c);
-        dboxes & b = boxes.at(ml).at(rl).at(c);
+        dboxes & b = boxes.AT(ml).AT(rl).AT(c);
         
         // 1. Remove all other interiors from this owned region
         for (int cc=0; cc<h.components(rl); ++cc) {
           if (cc != c) {
-            b.owned -= boxes.at(ml).at(rl).at(cc).interior;
+            b.owned -= boxes.AT(ml).AT(rl).AT(cc).interior;
           }
         }
         b.owned.normalize();
         
         // 2. Make disjoint from all earlier owned regions
         for (int cc=0; cc<c; ++cc) {
-          b.owned -= boxes.at(ml).at(rl).at(cc).owned;
+          b.owned -= boxes.AT(ml).AT(rl).AT(cc).owned;
         }
         b.owned.normalize();
         
@@ -201,7 +201,7 @@ void dh::foreach_reflevel_component_mglevel (dh::boxesop op)
   for (int ml=0; ml<h.mglevels(); ++ml) {
     for (int rl=0; rl<h.reflevels(); ++rl) {
       for (int c=0; c<h.components(rl); ++c) {
-        dboxes & box = boxes.at(ml).at(rl).at(c);
+        dboxes & box = boxes.AT(ml).AT(rl).AT(c);
         (this->*op)(box, rl, c, ml); // evaluate member function
       }
     }
@@ -238,7 +238,7 @@ void dh::setup_sync_boxes (dh::dboxes & box,
 
   // Sync boxes
   for (int cc=0; cc<h.components(rl); ++cc) {
-    dboxes & box1 = boxes.at(ml).at(rl).at(cc);
+    dboxes & box1 = boxes.AT(ml).AT(rl).AT(cc);
     ibset ovlp;
     if (cc != c) {
       // Do not sync from the same component
@@ -247,8 +247,8 @@ void dh::setup_sync_boxes (dh::dboxes & box,
       ovlp.normalize();
     }
     for (ibset::const_iterator b=ovlp.begin(); b!=ovlp.end(); ++b) {
-      box .recv_sync.at(cc).push_back(*b);
-      box1.send_sync.at(c).push_back(*b);
+      box .recv_sync.AT(cc).push_back(*b);
+      box1.send_sync.AT(c).push_back(*b);
     }
   }
 }
@@ -261,7 +261,7 @@ void dh::setup_multigrid_boxes (dh::dboxes & box,
 
   // Multigrid boxes
   if (ml>0) {
-    dboxes & bbox = boxes.at(ml-1).at(rl).at(c);
+    dboxes & bbox = boxes.AT(ml-1).AT(rl).AT(c);
     const ibbox intrf = bbox.interior;
     const ibbox extrf = bbox.exterior;
     // Restriction (interior)
@@ -303,7 +303,7 @@ void dh::setup_refinement_prolongation_boxes (dh::dboxes & box,
   // Refinement boxes
   if (rl<h.reflevels()-1) {
     for (int cc=0; cc<h.components(rl+1); ++cc) {
-      dboxes & box1 = boxes.at(ml).at(rl+1).at(cc);
+      dboxes & box1 = boxes.AT(ml).AT(rl+1).AT(cc);
       const ibbox intrf = box1.interior;
       // Prolongation (interior)
       // TODO: prefer boxes from the same processor
@@ -331,8 +331,8 @@ void dh::setup_refinement_prolongation_boxes (dh::dboxes & box,
           const ibbox send = recv.expanded_for(extr);
           assert (! send.empty());
           assert (send.is_contained_in(extr));
-          box1.recv_ref_coarse.at(c).push_back(recv);
-          box. send_ref_fine  .at(cc).push_back(send);
+          box1.recv_ref_coarse.AT(c).push_back(recv);
+          box. send_ref_fine  .AT(cc).push_back(send);
         }
       }
     } // for cc
@@ -347,7 +347,7 @@ void dh::setup_refinement_boundary_prolongation_boxes (dh::dboxes & box,
   // Refinement boxes
   if (rl<h.reflevels()-1) {
     for (int cc=0; cc<h.components(rl+1); ++cc) {
-      dboxes & box1 = boxes.at(ml).at(rl+1).at(cc);
+      dboxes & box1 = boxes.AT(ml).AT(rl+1).AT(cc);
       const ibbox & intrf = box1.interior;
       const ibbox & extrf = box1.exterior;
       const ibset & bndsf = box1.boundaries;
@@ -414,8 +414,8 @@ void dh::setup_refinement_boundary_prolongation_boxes (dh::dboxes & box,
             assert (! send.empty());
             assert (send.is_contained_in(extr));
             assert (send.is_contained_in(extr.expand(-pss,-pss)));
-            box1.recv_ref_bnd_coarse.at(c).push_back(recv);
-            box .send_ref_bnd_fine  .at(cc).push_back(send);
+            box1.recv_ref_bnd_coarse.AT(c).push_back(recv);
+            box .send_ref_bnd_fine  .AT(cc).push_back(send);
           }
         }
       }
@@ -434,7 +434,7 @@ void dh::setup_refinement_restriction_boxes (dh::dboxes & box,
   // Refinement boxes
   if (rl<h.reflevels()-1) {
     for (int cc=0; cc<h.components(rl+1); ++cc) {
-      dboxes & box1 = boxes.at(ml).at(rl+1).at(cc);
+      dboxes & box1 = boxes.AT(ml).AT(rl+1).AT(cc);
       const ibbox intrf = box1.interior;
       // Restriction (interior)
       {
@@ -471,7 +471,7 @@ void dh::setup_refinement_restriction_boxes (dh::dboxes & box,
           // remove what is sent during boundary prolongation
           const int pss = prolongation_stencil_size();
           for (int ccc=0; ccc<h.components(rl); ++ccc) {
-            const dh::dboxes& box2 = boxes.at(ml).at(rl).at(ccc);
+            const dh::dboxes& box2 = boxes.AT(ml).AT(rl).AT(ccc);
             for (iblistvect::const_iterator slvi =
                    box2.send_ref_bnd_fine.begin();
                  slvi != box2.send_ref_bnd_fine.end(); ++ slvi)
@@ -495,8 +495,8 @@ void dh::setup_refinement_restriction_boxes (dh::dboxes & box,
           assert (! recv.empty());
           const ibbox & send = recv.expanded_for(intrf);
           assert (! send.empty());
-          box1.send_ref_coarse.at(c).push_back(send);
-          box .recv_ref_fine  .at(cc).push_back(recv);
+          box1.send_ref_coarse.AT(c).push_back(send);
+          box .recv_ref_fine  .AT(cc).push_back(recv);
         }
       }
             
@@ -572,8 +572,8 @@ void dh::check_bboxes (dh::dboxes & box,
 #endif
     
     // Check that a component does not sync from itself
-    assert (box.recv_sync.at(c).empty());
-    assert (box.send_sync.at(c).empty());
+    assert (box.recv_sync.AT(c).empty());
+    assert (box.send_sync.AT(c).empty());
   }
 
   // Assert that the interior is received exactly once during
@@ -682,7 +682,7 @@ void dh::check_bboxes (dh::dboxes & box,
     ibset combined_extr;
     for (int cc=0; cc<h.components(rl); ++cc) {
       
-      dboxes const & box1 = boxes.at(ml).at(rl).at(cc);
+      dboxes const & box1 = boxes.AT(ml).AT(rl).AT(cc);
       
       ibbox const & intr = box1.interior;
       ibset const & ownd = box1.owned;
@@ -804,7 +804,7 @@ void dh::check_bboxes (dh::dboxes & box,
       // Points that are filled by prolongation from level rl-1
       ibset recvs;
       for (int cc=0; cc<h.components(rl); ++cc) {
-        dboxes & box1 = boxes.at(ml).at(rl).at(cc);
+        dboxes & box1 = boxes.AT(ml).AT(rl).AT(cc);
         for (iblistvect::const_iterator rlvi = box1.recv_ref_bnd_coarse.begin();
              rlvi != box1.recv_ref_bnd_coarse.end(); ++ rlvi)
         {
@@ -830,7 +830,7 @@ void dh::check_bboxes (dh::dboxes & box,
       // Points that are used for prolongation to level rl+1
       ibset sends;
       for (int cc=0; cc<h.components(rl); ++cc) {
-        dboxes & box1 = boxes.at(ml).at(rl).at(cc);
+        dboxes & box1 = boxes.AT(ml).AT(rl).AT(cc);
         for (iblistvect::const_iterator slvi = box1.send_ref_bnd_fine.begin();
              slvi != box1.send_ref_bnd_fine.end(); ++ slvi)
         {
@@ -865,13 +865,13 @@ void dh::calculate_bases ()
   // Calculate bases
   bases.resize(h.mglevels());
   for (int ml=0; ml<h.mglevels(); ++ml) {
-    bases.at(ml).resize(h.reflevels());
+    bases.AT(ml).resize(h.reflevels());
     for (int rl=0; rl<h.reflevels(); ++rl) {
-      dbases & base = bases.at(ml).at(rl);
+      dbases & base = bases.AT(ml).AT(rl);
       base.exterior = ibbox();
       base.interior = ibbox();
       for (int c=0; c<h.components(rl); ++c) {
-        dboxes & box = boxes.at(ml).at(rl).at(c);
+        dboxes & box = boxes.AT(ml).AT(rl).AT(c);
         base.exterior = (base.exterior.expanded_containing(box.exterior));
         base.interior = (base.interior.expanded_containing(box.interior));
       }
@@ -943,7 +943,7 @@ void dh::output_bases ()
   for (int ml=0; ml<h.mglevels(); ++ml) {
     for (int rl=0; rl<h.reflevels(); ++rl) {
       if (h.components(rl)>0) {
-        dbases & base = bases.at(ml).at(rl);
+        dbases & base = bases.AT(ml).AT(rl);
         cout << endl;
         cout << "dh bases:" << endl;
         cout << "ml=" << ml << " rl=" << rl << endl;
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index 207cc2624..7180c2ce9 100644
--- a/Carpet/CarpetLib/src/gdata.cc
+++ b/Carpet/CarpetLib/src/gdata.cc
@@ -128,15 +128,15 @@ void gdata::copy_from (comm_state& state,
     if (proc() != src->proc()) {
       // if this is a destination processor: advance its recv buffer size
       vector<comm_state::procbufdesc>& procbufs =
-        state.typebufs.at(c_datatype()).procbufs;
+        state.typebufs.AT(c_datatype()).procbufs;
       if (proc() == dist::rank()) {
-        procbufs.at(src->proc()).recvbufsize += box.size();
-        state.typebufs.at(c_datatype()).in_use = true;
+        procbufs.AT(src->proc()).recvbufsize += box.size();
+        state.typebufs.AT(c_datatype()).in_use = true;
       }
       // if this is a source processor: increment its send buffer size
       if (src->proc() == dist::rank()) {
-        procbufs.at(proc()).sendbufsize += box.size();
-        state.typebufs.at(c_datatype()).in_use = true;
+        procbufs.AT(proc()).sendbufsize += box.size();
+        state.typebufs.AT(c_datatype()).in_use = true;
       }
     }
     break;
@@ -157,7 +157,7 @@ void gdata::copy_from (comm_state& state,
     // if this is a destination processor and data has already been received
     // from the source processor: copy it from the recv buffer
     if (proc() == dist::rank() and
-        state.recvbuffers_ready.at(dist::size()*c_datatype() + src->proc())) {
+        state.recvbuffers_ready.AT(dist::size()*c_datatype() + src->proc())) {
       copy_from_recvbuffer (state, src, box);
     }
     break;
@@ -258,7 +258,7 @@ void gdata::copy_from_wait (comm_state& state,
     ivect items = (box.upper() - box.lower()) / box.stride() + 1;
     ivect offs  = (box.lower() - ext.lower()) / ext.stride();
     const char* recv_buffer = (const char*) b->pointer();
-    int& datatypesize = state.typebufs.at(c_datatype()).datatypesize;
+    int& datatypesize = state.typebufs.AT(c_datatype()).datatypesize;
 
     for (int k = 0; k < items[2]; k++) {
       for (int j = 0; j < items[1]; j++) {
@@ -293,9 +293,9 @@ void gdata::copy_into_sendbuffer (comm_state& state,
   } else {
     // copy to remote processor
     assert (src->_has_storage);
-    int datatypesize = state.typebufs.at(c_datatype()).datatypesize;
+    int const datatypesize = state.typebufs.AT(c_datatype()).datatypesize;
     comm_state::procbufdesc& procbuf =
-      state.typebufs.at(c_datatype()).procbufs.at(proc());
+      state.typebufs.AT(c_datatype()).procbufs.AT(proc());
     assert (procbuf.sendbuf - procbuf.sendbufbase <=
             ((int)procbuf.sendbufsize - box.size()) * datatypesize);
     int const fillstate = procbuf.sendbuf + (int)box.size()*datatypesize -
@@ -324,7 +324,7 @@ void gdata::copy_into_sendbuffer (comm_state& state,
       if (fillstate == (int)procbuf.sendbufsize * datatypesize) {
         wtime_commstate_isend.start();
         MPI_Isend (procbuf.sendbufbase, procbuf.sendbufsize,
-                   state.typebufs.at(c_datatype()).mpi_datatype,
+                   state.typebufs.AT(c_datatype()).mpi_datatype,
                    proc(), c_datatype(), dist::comm(),
                    &state.srequests.at(dist::size()*c_datatype() + proc()));
         wtime_commstate_isend.stop(procbuf.sendbufsize * datatypesize);
@@ -339,9 +339,9 @@ void gdata::copy_into_sendbuffer (comm_state& state,
 void gdata::copy_from_recvbuffer (comm_state& state,
                                   const gdata* src, const ibbox& box)
 {
-  int& datatypesize = state.typebufs.at(c_datatype()).datatypesize;
+  int& datatypesize = state.typebufs.AT(c_datatype()).datatypesize;
   comm_state::procbufdesc& procbuf =
-    state.typebufs.at(c_datatype()).procbufs.at(src->proc());
+    state.typebufs.AT(c_datatype()).procbufs.AT(src->proc());
   assert (procbuf.recvbuf - procbuf.recvbufbase <=
           ((int)procbuf.recvbufsize-box.size()) * datatypesize);
 
@@ -385,12 +385,12 @@ void gdata
   assert (all((box.lower()-extent().lower())%box.stride() == 0));
   assert (srcs.size() == times.size() and srcs.size()>0);
   for (int t=0; t<(int)srcs.size(); ++t) {
-    assert (srcs.at(t)->has_storage());
-    assert (all(box.lower()>=srcs.at(t)->extent().lower()));
-    assert (all(box.upper()<=srcs.at(t)->extent().upper()));
+    assert (srcs.AT(t)->has_storage());
+    assert (all(box.lower()>=srcs.AT(t)->extent().lower()));
+    assert (all(box.upper()<=srcs.AT(t)->extent().upper()));
   }
   assert (not box.empty());
-  const gdata* src = srcs.at(0);
+  const gdata* src = srcs.AT(0);
   if (dist::rank() != proc() and dist::rank() != src->proc()) return;
 
   switch (state.thestate) {
@@ -413,15 +413,15 @@ void gdata
     if (proc() != src->proc()) {
       // if this is a destination processor: increment its recv buffer size
       vector<comm_state::procbufdesc>& procbufs =
-        state.typebufs.at(c_datatype()).procbufs;
+        state.typebufs.AT(c_datatype()).procbufs;
       if (proc() == dist::rank()) {
-        procbufs.at(src->proc()).recvbufsize += box.size();
-        state.typebufs.at(c_datatype()).in_use = true;
+        procbufs.AT(src->proc()).recvbufsize += box.size();
+        state.typebufs.AT(c_datatype()).in_use = true;
       }
       // if this is a source processor: increment its send buffer size
       if (src->proc() == dist::rank()) {
-        procbufs.at(proc()).sendbufsize += box.size();
-        state.typebufs.at(c_datatype()).in_use = true;
+        procbufs.AT(proc()).sendbufsize += box.size();
+        state.typebufs.AT(c_datatype()).in_use = true;
       }
     }
     break;
@@ -443,7 +443,7 @@ void gdata
     // from the source processor: copy it from the recv buffer
     // (the processor-local case is not handled here)
     if (proc() == dist::rank() and
-        state.recvbuffers_ready.at(dist::size()*c_datatype() + src->proc())) {
+        state.recvbuffers_ready.AT(dist::size()*c_datatype() + src->proc())) {
       copy_from_recvbuffer(state, src, box);
     }
     break;
@@ -462,7 +462,7 @@ void gdata
                          const int order_space,
                          const int order_time)
 {
-  const gdata* src = srcs.at(0);
+  const gdata* src = srcs.AT(0);
   if (dist::rank() == proc()) {
     // interpolate from other processor
 
@@ -514,17 +514,17 @@ void gdata
 {
   DECLARE_CCTK_PARAMETERS;
   
-  if (proc() == srcs.at(0)->proc()) {
+  if (proc() == srcs.AT(0)->proc()) {
     // interpolate on same processor
     interpolate_from_innerloop (srcs, times, box, time,
                                 order_space, order_time);
   } else {
     // interpolate to remote processor
-    const gdata* src = srcs.at(0);
+    const gdata* src = srcs.AT(0);
     assert (src->_has_storage);
-    int& datatypesize = state.typebufs.at(c_datatype()).datatypesize;
+    int& datatypesize = state.typebufs.AT(c_datatype()).datatypesize;
     comm_state::procbufdesc& procbuf =
-      state.typebufs.at(c_datatype()).procbufs.at(proc());
+      state.typebufs.AT(c_datatype()).procbufs.AT(proc());
     assert (procbuf.sendbuf - procbuf.sendbufbase <=
             ((int)procbuf.sendbufsize - box.size()) * datatypesize);
     assert (src->has_storage());
diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc
index 26fa00f38..50bb0f661 100644
--- a/Carpet/CarpetLib/src/gf.cc
+++ b/Carpet/CarpetLib/src/gf.cc
@@ -58,7 +58,7 @@ const data<T>* gf<T>::operator() (int tl, int rl, int c, int ml) const
   assert (c>=0 and c<h.components(rl));
   assert (ml>=0 and ml<h.mglevels());
   assert (tl>=0 and tl<timelevels(ml, rl));
-  return (const data<T>*)storage.at(ml).at(rl).at(c).at(tl);
+  return (const data<T>*)storage.AT(ml).AT(rl).AT(c).AT(tl);
 }
 
 template<typename T>
@@ -68,7 +68,7 @@ data<T>* gf<T>::operator() (int tl, int rl, int c, int ml)
   assert (c>=0 and c<h.components(rl));
   assert (ml>=0 and ml<h.mglevels());
   assert (tl>=0 and tl<timelevels(ml, rl));
-  return (data<T>*)storage.at(ml).at(rl).at(c).at(tl);
+  return (data<T>*)storage.AT(ml).AT(rl).AT(c).AT(tl);
 }
 
 
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index e81a3811f..895cf941d 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -38,7 +38,7 @@ ggf::ggf (const int varindex_, const operator_type transport_operator_,
   
   timelevels_.resize(d.h.mglevels());
   for (int ml=0; ml<d.h.mglevels(); ++ml) {
-    timelevels_.at(ml).resize(d.h.reflevels(), 0);
+    timelevels_.AT(ml).resize(d.h.reflevels(), 0);
   }
   
   d.add(this);
@@ -60,33 +60,33 @@ bool ggf::operator== (const ggf& f) const {
 void ggf::set_timelevels (const int ml, const int rl, const int new_timelevels)
 {
   assert (ml>=0 and ml<(int)storage.size());
-  assert (rl>=0 and rl<(int)storage.at(ml).size());
+  assert (rl>=0 and rl<(int)storage.AT(ml).size());
   
   assert (new_timelevels >= 0);
   
   if (new_timelevels < timelevels(ml,rl)) {
     
-    for (int c=0; c<(int)storage.at(ml).at(rl).size(); ++c) {
+    for (int c=0; c<(int)storage.AT(ml).AT(rl).size(); ++c) {
       for (int tl=new_timelevels; tl<timelevels(ml,rl); ++tl) {
-        delete storage.at(ml).at(rl).at(c).at(tl);
+        delete storage.AT(ml).AT(rl).AT(c).AT(tl);
       }
-      storage.at(ml).at(rl).at(c).resize (new_timelevels);
+      storage.AT(ml).AT(rl).AT(c).resize (new_timelevels);
     } // for c
     
   } else if (new_timelevels > timelevels(ml,rl)) {
     
-    for (int c=0; c<(int)storage.at(ml).at(rl).size(); ++c) {
-      storage.at(ml).at(rl).at(c).resize (new_timelevels);
+    for (int c=0; c<(int)storage.AT(ml).AT(rl).size(); ++c) {
+      storage.AT(ml).AT(rl).AT(c).resize (new_timelevels);
       for (int tl=timelevels(ml,rl); tl<new_timelevels; ++tl) {
-        storage.at(ml).at(rl).at(c).at(tl) = typed_data(tl,rl,c,ml);
-        storage.at(ml).at(rl).at(c).at(tl)->allocate
-          (d.boxes.at(ml).at(rl).at(c).exterior, h.processor(rl,c));
+        storage.AT(ml).AT(rl).AT(c).AT(tl) = typed_data(tl,rl,c,ml);
+        storage.AT(ml).AT(rl).AT(c).AT(tl)->allocate
+          (d.boxes.AT(ml).AT(rl).AT(c).exterior, h.processor(rl,c));
       } // for tl
     } // for c
     
   }
   
-  timelevels_.at(ml).at(rl) = new_timelevels;
+  timelevels_.AT(ml).AT(rl) = new_timelevels;
 }
 
 
@@ -95,14 +95,14 @@ void ggf::recompose_crop ()
 {
   // Free storage that will not be needed
   for (int ml=0; ml<h.mglevels(); ++ml) {
-    for (int rl=h.reflevels(); rl<(int)storage.at(ml).size(); ++rl) {
-      for (int c=0; c<(int)storage.at(ml).at(rl).size(); ++c) {
-        for (int tl=0; tl<(int)storage.at(ml).at(rl).at(c).size(); ++tl) {
-          delete storage.at(ml).at(rl).at(c).at(tl);
+    for (int rl=h.reflevels(); rl<(int)storage.AT(ml).size(); ++rl) {
+      for (int c=0; c<(int)storage.AT(ml).AT(rl).size(); ++c) {
+        for (int tl=0; tl<(int)storage.AT(ml).AT(rl).AT(c).size(); ++tl) {
+          delete storage.AT(ml).AT(rl).AT(c).AT(tl);
         } // for tl
       } // for c
     } // for rl
-    storage.at(ml).resize(h.reflevels());
+    storage.AT(ml).resize(h.reflevels());
   } // for ml
 }
 
@@ -111,26 +111,26 @@ void ggf::recompose_allocate (const int rl)
   // Retain storage that might be needed
   oldstorage.resize(storage.size());
   for (int ml=0; ml<(int)storage.size(); ++ml) {
-    oldstorage.at(ml).resize(storage.at(ml).size());
-    oldstorage.at(ml).at(rl) = storage.at(ml).at(rl);
-    storage.at(ml).at(rl).resize(0);
+    oldstorage.AT(ml).resize(storage.AT(ml).size());
+    oldstorage.AT(ml).AT(rl) = storage.AT(ml).AT(rl);
+    storage.AT(ml).AT(rl).resize(0);
   }
   
   for (int ml=0; ml<d.h.mglevels(); ++ml) {
-    timelevels_.at(ml).resize(d.h.reflevels(), timelevels_.at(ml).at(0));
+    timelevels_.AT(ml).resize(d.h.reflevels(), timelevels_.AT(ml).AT(0));
   }
   
   // Resize structure and allocate storage
   storage.resize(h.mglevels());
   for (int ml=0; ml<h.mglevels(); ++ml) {
-    storage.at(ml).resize(h.reflevels());
-    storage.at(ml).at(rl).resize(h.components(rl));
+    storage.AT(ml).resize(h.reflevels());
+    storage.AT(ml).AT(rl).resize(h.components(rl));
     for (int c=0; c<h.components(rl); ++c) {
-      storage.at(ml).at(rl).at(c).resize(timelevels(ml,rl));
+      storage.AT(ml).AT(rl).AT(c).resize(timelevels(ml,rl));
       for (int tl=0; tl<timelevels(ml,rl); ++tl) {
-        storage.at(ml).at(rl).at(c).at(tl) = typed_data(tl,rl,c,ml);
-        storage.at(ml).at(rl).at(c).at(tl)->allocate
-          (d.boxes.at(ml).at(rl).at(c).exterior, h.processor(rl,c));
+        storage.AT(ml).AT(rl).AT(c).AT(tl) = typed_data(tl,rl,c,ml);
+        storage.AT(ml).AT(rl).AT(c).AT(tl)->allocate
+          (d.boxes.AT(ml).AT(rl).AT(c).exterior, h.processor(rl,c));
       } // for tl
     } // for c
   } // for ml
@@ -147,20 +147,20 @@ void ggf::recompose_fill (comm_state& state, const int rl,
         // Find out which regions need to be prolongated
         // (Copy the exterior because some variables are not prolongated)
         // TODO: do this once in the dh instead of for each variable here
-        ibset work (d.boxes.at(ml).at(rl).at(c).exterior);
+        ibset work (d.boxes.AT(ml).AT(rl).AT(c).exterior);
         
         // Copy from old storage, if possible
         // TODO: copy only from interior regions?
-        if (rl<(int)oldstorage.at(ml).size()) {
-          for (int cc=0; cc<(int)oldstorage.at(ml).at(rl).size(); ++cc) {
+        if (rl<(int)oldstorage.AT(ml).size()) {
+          for (int cc=0; cc<(int)oldstorage.AT(ml).AT(rl).size(); ++cc) {
             // TODO: prefer same processor, etc., see dh.cc
             ibset ovlp
-              = work & oldstorage.at(ml).at(rl).at(cc).at(tl)->extent();
+              = work & oldstorage.AT(ml).AT(rl).AT(cc).AT(tl)->extent();
             ovlp.normalize();
             work -= ovlp;
             for (ibset::const_iterator r=ovlp.begin(); r!=ovlp.end(); ++r) {
-              storage.at(ml).at(rl).at(c).at(tl)->copy_from
-                (state, oldstorage.at(ml).at(rl).at(cc).at(tl), *r);
+              storage.AT(ml).AT(rl).AT(c).AT(tl)->copy_from
+                (state, oldstorage.AT(ml).AT(rl).AT(cc).AT(tl), *r);
             }
           } // for cc
         } // if rl
@@ -178,21 +178,21 @@ void ggf::recompose_fill (comm_state& state, const int rl,
               vector<int> tls(numtl);
               vector<CCTK_REAL> times(numtl);
               for (int i=0; i<numtl; ++i) {
-                tls.at(i) = i;
-                times.at(i) = t.time(tls.at(i),rl-1,ml);
+                tls.AT(i) = i;
+                times.AT(i) = t.time(tls.AT(i),rl-1,ml);
               }
-              for (int cc=0; cc<(int)storage.at(ml).at(rl-1).size(); ++cc) {
+              for (int cc=0; cc<(int)storage.AT(ml).AT(rl-1).size(); ++cc) {
                 vector<const gdata*> gsrcs(numtl);
                 for (int i=0; i<numtl; ++i) {
-                  gsrcs.at(i) = storage.at(ml).at(rl-1).at(cc).at(tls.at(i));
-                  assert (gsrcs.at(i)->extent() == gsrcs.at(0)->extent());
+                  gsrcs.AT(i) = storage.AT(ml).AT(rl-1).AT(cc).AT(tls.AT(i));
+                  assert (gsrcs.AT(i)->extent() == gsrcs.AT(0)->extent());
                 }
                 const CCTK_REAL time = t.time(tl,rl,ml);
                 
                 // TODO: choose larger regions first
                 // TODO: prefer regions from the same processor
                 const iblist& list
-                  = d.boxes.at(ml).at(rl).at(c).recv_ref_coarse.at(cc);
+                  = d.boxes.AT(ml).AT(rl).AT(c).recv_ref_coarse.AT(cc);
                 for (iblist::const_iterator iter=list.begin();
                      iter!=list.end(); ++iter)
                 {
@@ -202,7 +202,7 @@ void ggf::recompose_fill (comm_state& state, const int rl,
                   for (ibset::const_iterator r=ovlp.begin();
                        r!=ovlp.end(); ++r)
                   {
-                    storage.at(ml).at(rl).at(c).at(tl)->interpolate_from
+                    storage.AT(ml).AT(rl).AT(c).AT(tl)->interpolate_from
                       (state, gsrcs, times, *r, time, pos, pot);
                   } // for r
                 } // for iter
@@ -225,12 +225,12 @@ void ggf::recompose_free (const int rl)
 {
   // Delete old storage
   for (int ml=0; ml<(int)oldstorage.size(); ++ml) {
-    for (int c=0; c<(int)oldstorage.at(ml).at(rl).size(); ++c) {
+    for (int c=0; c<(int)oldstorage.AT(ml).AT(rl).size(); ++c) {
       for (int tl=0; tl<timelevels(ml,rl); ++tl) {
-        delete oldstorage.at(ml).at(rl).at(c).at(tl);
+        delete oldstorage.AT(ml).AT(rl).AT(c).AT(tl);
       } // for tl
     } // for c
-    oldstorage.at(ml).at(rl).resize(0);
+    oldstorage.AT(ml).AT(rl).resize(0);
   } // for ml
 }
 
@@ -295,9 +295,9 @@ void ggf::flip (int rl, int c, int ml) {
     const int tl1 =                  tl;
     const int tl2 = timelevels(ml,rl)-1 - tl;
     assert (tl1 < tl2);
-    gdata* tmpdata = storage.at(ml).at(rl).at(c).at(tl1);
-    storage.at(ml).at(rl).at(c).at(tl1) = storage.at(ml).at(rl).at(c).at(tl2);
-    storage.at(ml).at(rl).at(c).at(tl2) = tmpdata;
+    gdata* tmpdata = storage.AT(ml).AT(rl).AT(c).AT(tl1);
+    storage.AT(ml).AT(rl).AT(c).AT(tl1) = storage.AT(ml).AT(rl).AT(c).AT(tl2);
+    storage.AT(ml).AT(rl).AT(c).AT(tl2) = tmpdata;
   }
 }
 
@@ -316,13 +316,13 @@ void ggf::copycat (comm_state& state,
   assert (ml1>=0 and ml1<h.mglevels());
   assert (tl1>=0 and tl1<timelevels(ml1,rl1));
   assert (rl2>=0 and rl2<h.reflevels());
-  const int c2=c1;
+  int const c2= c1;
   assert (ml2<h.mglevels());
   assert (tl2>=0 and tl2<timelevels(ml2,rl2));
-  const ibbox recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_box;
+  ibbox const & recv = d.boxes.AT(ml1).AT(rl1).AT(c1).*recv_box;
   // copy the content
-  gdata* const dst = storage.at(ml1).at(rl1).at(c1).at(tl1);
-  gdata* const src = storage.at(ml2).at(rl2).at(c2).at(tl2);
+  gdata * const dst = storage.AT(ml1).AT(rl1).AT(c1).AT(tl1);
+  gdata const * const src = storage.AT(ml2).AT(rl2).AT(c2).AT(tl2);
   dst->copy_from(state, src, recv);
 }
 
@@ -338,15 +338,15 @@ void ggf::copycat (comm_state& state,
   assert (tl1>=0 and tl1<timelevels(ml1,rl1));
   assert (           ml2<h.mglevels());
   assert (rl2>=0 and rl2<h.reflevels());
-  const int c2=c1;
+  int const c2 = c1;
   assert (tl2>=0 and tl2<timelevels(ml2,rl2));
-  const iblist recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_list;
+  iblist const & recv = d.boxes.AT(ml1).AT(rl1).AT(c1).*recv_list;
   // walk all boxes
   for (iblist::const_iterator r=recv.begin(); r!=recv.end(); ++r) {
     // (use the send boxes for communication)
     // copy the content
-    gdata* const dst = storage.at(ml1).at(rl1).at(c1).at(tl1);
-    gdata* const src = storage.at(ml2).at(rl2).at(c2).at(tl2);
+    gdata * const dst = storage.AT(ml1).AT(rl1).AT(c1).AT(tl1);
+    gdata const * const src = storage.AT(ml2).AT(rl2).AT(c2).AT(tl2);
     dst->copy_from(state, src, *r);
   }
 }
@@ -390,22 +390,22 @@ void ggf::intercat (comm_state& state,
   assert (ml1>=0 and ml1<h.mglevels());
   assert (tl1>=0 and tl1<timelevels(ml1,rl1));
   assert (rl2>=0 and rl2<h.reflevels());
-  const int c2=c1;
+  int const c2 = c1;
   assert (ml2>=0 and ml2<h.mglevels());
   for (int i=0; i<(int)tl2s.size(); ++i) {
-    assert (tl2s.at(i)>=0 and tl2s.at(i)<timelevels(ml2,rl2));
+    assert (tl2s.AT(i)>=0 and tl2s.AT(i)<timelevels(ml2,rl2));
   }
   
   vector<const gdata*> gsrcs(tl2s.size());
   vector<CCTK_REAL> times(tl2s.size());
   for (int i=0; i<(int)gsrcs.size(); ++i) {
-    gsrcs.at(i) = storage.at(ml2).at(rl2).at(c2).at(tl2s.at(i));
-    times.at(i) = t.time(tl2s.at(i),rl2,ml2);
+    gsrcs.AT(i) = storage.AT(ml2).AT(rl2).AT(c2).AT(tl2s.AT(i));
+    times.AT(i) = t.time(tl2s.AT(i),rl2,ml2);
   }
   
-  const ibbox recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_list;
+  ibbox const & recv = d.boxes.AT(ml1).AT(rl1).AT(c1).*recv_list;
   // interpolate the content
-  storage.at(ml1).at(rl1).at(c1).at(tl1)->interpolate_from
+  storage.AT(ml1).AT(rl1).AT(c1).AT(tl1)->interpolate_from
     (state, gsrcs, times, recv, time,
      d.prolongation_order_space, prolongation_order_time);
 }
@@ -422,26 +422,26 @@ void ggf::intercat (comm_state& state,
   assert (ml1>=0 and ml1<h.mglevels());
   assert (tl1>=0 and tl1<timelevels(ml1,rl1));
   assert (rl2>=0 and rl2<h.reflevels());
-  const int c2=c1;
+  int const c2 = c1;
   assert (ml2>=0 and ml2<h.mglevels());
   for (int i=0; i<(int)tl2s.size(); ++i) {
-    assert (tl2s.at(i)>=0 and tl2s.at(i)<timelevels(ml2,rl2));
+    assert (tl2s.AT(i)>=0 and tl2s.AT(i)<timelevels(ml2,rl2));
   }
   
   vector<const gdata*> gsrcs(tl2s.size());
   vector<CCTK_REAL> times(tl2s.size());
   for (int i=0; i<(int)gsrcs.size(); ++i) {
-    gsrcs.at(i) = storage.at(ml2).at(rl2).at(c2).at(tl2s.at(i));
-    times.at(i) = t.time(tl2s.at(i),rl2,ml2);
+    gsrcs.AT(i) = storage.AT(ml2).AT(rl2).AT(c2).AT(tl2s.AT(i));
+    times.AT(i) = t.time(tl2s.AT(i),rl2,ml2);
   }
   
-  const iblist recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_list;
+  iblist const & recv = d.boxes.AT(ml1).AT(rl1).AT(c1).*recv_list;
   // walk all boxes
   for (iblist::const_iterator r=recv.begin(); r!=recv.end(); ++r)
   {
     // (use the send boxes for communication)
     // interpolate the content
-    storage.at(ml1).at(rl1).at(c1).at(tl1)->interpolate_from
+    storage.AT(ml1).AT(rl1).AT(c1).AT(tl1)->interpolate_from
       (state, gsrcs, times, *r, time,
        d.prolongation_order_space, prolongation_order_time);
   }
@@ -460,7 +460,7 @@ void ggf::intercat (comm_state& state,
   assert (tl1>=0 and tl1<timelevels(ml1,rl1));
   assert (rl2>=0 and rl2<h.reflevels());
   for (int i=0; i<(int)tl2s.size(); ++i) {
-    assert (tl2s.at(i)>=0 and tl2s.at(i)<timelevels(ml2,rl2));
+    assert (tl2s.AT(i)>=0 and tl2s.AT(i)<timelevels(ml2,rl2));
   }
   // walk all components
   for (int c2=0; c2<h.components(rl2); ++c2) {
@@ -469,20 +469,21 @@ void ggf::intercat (comm_state& state,
     vector<const gdata*> gsrcs(tl2s.size());
     vector<CCTK_REAL> times(tl2s.size());
     for (int i=0; i<(int)gsrcs.size(); ++i) {
-      gsrcs.at(i) = storage.at(ml2).at(rl2).at(c2).at(tl2s.at(i));
-      times.at(i) = t.time(tl2s.at(i),rl2,ml2);
+      gsrcs.AT(i) = storage.AT(ml2).AT(rl2).AT(c2).AT(tl2s.AT(i));
+      times.AT(i) = t.time(tl2s.AT(i),rl2,ml2);
     }
     
-    const iblist recv = (d.boxes.at(ml1).at(rl1).at(c1).*recv_listvect).at(c2);
+    iblist const &recv =
+      (d.boxes.AT(ml1).AT(rl1).AT(c1).*recv_listvect).AT(c2);
     // walk all boxes
     for (iblist::const_iterator r=recv.begin(); r!=recv.end(); ++r)
     {
       // (use the send boxes for communication)
       // interpolate the content
       int const pos = d.prolongation_order_space;
-      int const pot
-        = transport_operator == op_copy ? 0 : prolongation_order_time;
-      storage.at(ml1).at(rl1).at(c1).at(tl1)->interpolate_from
+      int const pot =
+        transport_operator == op_copy ? 0 : prolongation_order_time;
+      storage.AT(ml1).AT(rl1).AT(c1).AT(tl1)->interpolate_from
       	(state, gsrcs, times, *r, time, pos, pot);
     }
   }
@@ -529,11 +530,11 @@ void ggf::ref_bnd_prolongate (comm_state& state,
     }
     assert (timelevels(ml,rl) >= prolongation_order_time+1);
     tl2s.resize(prolongation_order_time+1);
-    for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = i;
+    for (int i=0; i<=prolongation_order_time; ++i) tl2s.AT(i) = i;
   } else {
     assert (timelevels(ml,rl) >= 1);
     tl2s.resize(1);
-    tl2s.at(0) = 0;
+    tl2s.AT(0) = 0;
   }
   intercat (state,
             tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_bnd_coarse,
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 9e47dfcff..3163983af 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -82,7 +82,7 @@ public:
   // Querying
   int timelevels (int const ml, int const rl) const
   {
-    return timelevels_.at(ml).at(rl);
+    return timelevels_.AT(ml).AT(rl);
   }
   
   
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index bb9f96e7d..159aa1e3d 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -28,10 +28,10 @@ gh::gh (const vector<ivect> & reffacts_, const centering refcent_,
   assert (reffacts.size() >= 1);
   assert (all (reffacts.front() == 1));
   for (size_t n = 1; n < reffacts.size(); ++ n) {
-    assert (all (reffacts.at(n) >= reffacts.at(n-1)));
-    assert (all (reffacts.at(n) % reffacts.at(n-1) == 0));
+    assert (all (reffacts.AT(n) >= reffacts.AT(n-1)));
+    assert (all (reffacts.AT(n) % reffacts.AT(n-1) == 0));
   }
-  assert (all (baseextent.stride() % reffacts.at(reffacts.size()-1) == 0));
+  assert (all (baseextent.stride() % reffacts.AT(reffacts.size()-1) == 0));
 }
 
 // Destructors
@@ -91,8 +91,8 @@ bool gh::recompose (const int rl,
     
     // Overwrite old with new grid hierarchy
     for (int ml=0; ml<mglevels(); ++ml) {
-      _oldregions.at(ml).resize (_regions.at(ml).size());
-      _oldregions.at(ml).at(rl) = _regions.at(ml).at(rl);
+      _oldregions.AT(ml).resize (_regions.AT(ml).size());
+      _oldregions.AT(ml).AT(rl) = _regions.AT(ml).AT(rl);
     }
     
   }
@@ -106,12 +106,12 @@ bool gh::level_did_change (const int rl) const
   if (_regions.size() != _oldregions.size()) return true;
   for (int ml=0; ml<mglevels(); ++ml) {
     assert (rl>=0 and rl<reflevels());
-    if (rl >= (int)_oldregions.at(ml).size()) return true;
-    if (_regions.at(ml).at(rl).size() != _oldregions.at(ml).at(rl).size()) {
+    if (rl >= (int)_oldregions.AT(ml).size()) return true;
+    if (_regions.AT(ml).AT(rl).size() != _oldregions.AT(ml).AT(rl).size()) {
       return true;
     }
     for (int c=0; c<components(rl); ++c) {
-      if (_regions.at(ml).at(rl).at(c) != _oldregions.at(ml).at(rl).at(c)) {
+      if (_regions.AT(ml).AT(rl).AT(c) != _oldregions.AT(ml).AT(rl).AT(c)) {
         return true;
       }
     } // for c
@@ -178,7 +178,7 @@ void gh::check_refinement_levels ()
   for (int ml=0; ml<mglevels(); ++ml) {
     for (int rl=1; rl<reflevels(); ++rl) {
       assert (all(extent(ml,rl-1,0).stride()
-                  == ((reffacts.at(rl) / reffacts.at(rl-1))
+                  == ((reffacts.AT(rl) / reffacts.AT(rl-1))
                       * extent(ml,rl,0).stride())));
       // Check contained-ness:
       // first take all coarse grids ...
@@ -228,10 +228,10 @@ void gh::calculate_base_extents_of_all_levels ()
 {
   _bases.resize(mglevels());
   for (int ml=0; ml<mglevels(); ++ml) {
-    _bases.at(ml).resize(reflevels());
+    _bases.AT(ml).resize(reflevels());
     for (int rl=0; rl<reflevels(); ++rl) {
-      _bases.at(ml).at(rl) = ibbox();
-      ibbox &bb = _bases.at(ml).at(rl);
+      _bases.AT(ml).AT(rl) = ibbox();
+      ibbox &bb = _bases.AT(ml).AT(rl);
       for (int c=0; c<components(rl); ++c) {
         bb = bb.expanded_containing(extent(ml,rl,c));
       }
@@ -301,7 +301,7 @@ void gh::do_output_bases (ostream& os) const
         os << endl;
         os << "gh bases:" << endl;
         os << "ml=" << ml << " rl=" << rl << endl;
-        os << "base=" << bases().at(ml).at(rl) << endl;
+        os << "base=" << bases().AT(ml).AT(rl) << endl;
       }
     }
   }
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index f97e15ef7..c48f3cd9c 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -92,22 +92,22 @@ public:
   
   ibbox extent (const int m, const int rl, const int c) const
   {
-    return _regions.at(m).at(rl).at(c).extent;
+    return _regions.AT(m).AT(rl).AT(c).extent;
   }
   
   b2vect outer_boundaries (const int rl, const int c) const
   {
-    return _regions.at(0).at(rl).at(c).outer_boundaries;
+    return _regions.AT(0).AT(rl).AT(c).outer_boundaries;
   }
 
   b2vect refinement_boundaries (const int rl, const int c) const
   {
-    return _regions.at(0).at(rl).at(c).refinement_boundaries;
+    return _regions.AT(0).AT(rl).AT(c).refinement_boundaries;
   }
 
   int processor (const int rl, const int c) const
   {
-    return _regions.at(0).at(rl).at(c).processor;
+    return _regions.AT(0).AT(rl).AT(c).processor;
   }
 
   int mglevels () const
@@ -118,12 +118,12 @@ public:
   int reflevels () const
   {
     if (mglevels() == 0) return 0;
-    return (int)_regions.at(0).size();
+    return (int)_regions.AT(0).size();
   }
   
   int components (const int rl) const
   {
-    return (int)_regions.at(0).at(rl).size();
+    return (int)_regions.AT(0).AT(rl).size();
   }
 
   bool is_local (const int rl, const int c) const
diff --git a/Carpet/CarpetLib/src/mem.cc b/Carpet/CarpetLib/src/mem.cc
index 0f4c3ac06..8776e642b 100644
--- a/Carpet/CarpetLib/src/mem.cc
+++ b/Carpet/CarpetLib/src/mem.cc
@@ -109,8 +109,8 @@ void mem<T>::
 register_client (size_t const vectorindex)
 {
   assert (vectorindex < vectorlength_);
-  assert (! clients_.at(vectorindex));
-  clients_.at(vectorindex) = true;
+  assert (! clients_.AT(vectorindex));
+  clients_.AT(vectorindex) = true;
   ++ num_clients_;
 }
 
@@ -119,8 +119,8 @@ void mem<T>::
 unregister_client (size_t const vectorindex)
 {
   assert (vectorindex < vectorlength_);
-  assert (clients_.at(vectorindex));
-  clients_.at(vectorindex) = false;
+  assert (clients_.AT(vectorindex));
+  clients_.AT(vectorindex) = false;
   assert (num_clients_ > 0);
   -- num_clients_;
 }
diff --git a/Carpet/CarpetLib/src/mem.hh b/Carpet/CarpetLib/src/mem.hh
index 5649bc250..7ccbd6d72 100644
--- a/Carpet/CarpetLib/src/mem.hh
+++ b/Carpet/CarpetLib/src/mem.hh
@@ -25,7 +25,7 @@ public:
   T * storage (size_t vectorindex) const
   {
     assert (vectorindex < vectorlength_);
-    assert (clients_.at(vectorindex));
+    assert (clients_.AT(vectorindex));
     return & storage_ [vectorindex * nelems_];
   }
   
diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc
index bf69d1b5c..e8d5bbb49 100644
--- a/Carpet/CarpetLib/src/th.cc
+++ b/Carpet/CarpetLib/src/th.cc
@@ -21,8 +21,8 @@ th::th (gh& h_, const vector<int> & reffacts_, const CCTK_REAL basedelta)
   assert (reffacts.size() >= 1);
   assert (reffacts.front() == 1);
   for (size_t n = 1; n < reffacts.size(); ++ n) {
-    assert (reffacts.at(n) >= reffacts.at(n-1));
-    assert (reffacts.at(n) % reffacts.at(n-1) == 0);
+    assert (reffacts.AT(n) >= reffacts.AT(n-1));
+    assert (reffacts.AT(n) % reffacts.AT(n-1) == 0);
   }
   h.add(this);
 }
@@ -40,21 +40,21 @@ void th::regrid ()
   times.resize(h.mglevels());
   deltas.resize(h.mglevels());
   for (int ml=0; ml<h.mglevels(); ++ml) {
-    const int old_reflevels = times.at(ml).size();
-    times.at(ml).resize(h.reflevels());
-    deltas.at(ml).resize(h.reflevels());
+    const int old_reflevels = times.AT(ml).size();
+    times.AT(ml).resize(h.reflevels());
+    deltas.AT(ml).resize(h.reflevels());
     for (int rl=0; rl<h.reflevels(); ++rl) {
       if (old_mglevels==0) {
-        times.at(ml).at(rl) = 0;
+        times.AT(ml).AT(rl) = 0;
       } else if (rl < old_reflevels) {
         // do nothing
       } else {
-        times.at(ml).at(rl) = times.at(ml).at(rl-1);
+        times.AT(ml).AT(rl) = times.AT(ml).AT(rl-1);
       }
       if (ml==0) {
-	deltas.at(ml).at(rl) = delta / reffacts.at(rl);
+	deltas.AT(ml).AT(rl) = delta / reffacts.AT(rl);
       } else {
-	deltas.at(ml).at(rl) = deltas.at(ml-1).at(rl) * h.mgfact;
+	deltas.AT(ml).AT(rl) = deltas.AT(ml-1).AT(rl) * h.mgfact;
       }
     }
   }
@@ -72,7 +72,7 @@ void th::output (ostream& os) const
     for (int rl=0; rl<h.reflevels(); ++rl) {
       os << sep
          << ml << ":" << rl << ":"
-	 << times.at(ml).at(rl) << "(" << deltas.at(ml).at(rl) << ")";
+	 << times.AT(ml).AT(rl) << "(" << deltas.AT(ml).AT(rl) << ")";
       sep = ",";
     }
   }
diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh
index 9772afaa5..b3dc315b0 100644
--- a/Carpet/CarpetLib/src/th.hh
+++ b/Carpet/CarpetLib/src/th.hh
@@ -54,14 +54,14 @@ public:
   {
     assert (rl>=0 and rl<h.reflevels());
     assert (ml>=0 and ml<h.mglevels());
-    return times.at(ml).at(rl);
+    return times.AT(ml).AT(rl);
   }
   
   void set_time (const int rl, const int ml, const CCTK_REAL t)
   {
     assert (rl>=0 and rl<h.reflevels());
     assert (ml>=0 and ml<h.mglevels());
-    times.at(ml).at(rl) = t;
+    times.AT(ml).AT(rl) = t;
   }
   
   void advance_time (const int rl, const int ml)
@@ -73,14 +73,14 @@ public:
   {
     assert (rl>=0 and rl<h.reflevels());
     assert (ml>=0 and ml<h.mglevels());
-    return deltas.at(ml).at(rl);
+    return deltas.AT(ml).AT(rl);
   }
   
   void set_delta (const int rl, const int ml, const CCTK_REAL dt)
   {
     assert (rl>=0 and rl<h.reflevels());
     assert (ml>=0 and ml<h.mglevels());
-    deltas.at(ml).at(rl) = dt;
+    deltas.AT(ml).AT(rl) = dt;
   }
   
   CCTK_REAL time (const int tl, const int rl, const int ml) const