diff --git a/include/ftk/mesh/aabb.hh b/include/ftk/mesh/aabb.hh
index 51477531..3ad7bbd2 100644
--- a/include/ftk/mesh/aabb.hh
+++ b/include/ftk/mesh/aabb.hh
@@ -15,25 +15,52 @@ inline T max3(T x, T y, T z) {
   return std::max(std::max(x, y), z);
 }
 
-template <typename T>
+template <int N, typename T>
 struct AABB {
   int id = 0;
-  T A[2] = {std::numeric_limits<T>::max(), std::numeric_limits<T>::max()}, 
-    B[2] = {-std::numeric_limits<T>::max(), -std::numeric_limits<T>::max()};
-  T C[2] = {0, 0}; // centroid
+  std::array<T, N> A, B, C;
+
+  // T A[2] = {std::numeric_limits<T>::max(), std::numeric_limits<T>::max()}, 
+  //   B[2] = {-std::numeric_limits<T>::max(), -std::numeric_limits<T>::max()};
+  // T C[2] = {0, 0}; // centroid
+
+  AABB() {
+    for (int i = 0; i < N; i ++) {
+      A[i] =  std::numeric_limits<T>::max();
+      B[i] = -std::numeric_limits<T>::max();
+      C[i] = T(0);
+    }
+  }
 
+  bool contains(const std::array<T, N> &X) const {
+    for (int i = 0; i < N; i ++) 
+      if (X[i] < A[i] || X[i] >= B[i])
+        return false;
+    return true;
+  }
+  
   bool contains(const T X[]) const {
-    return X[0] >= A[0] && X[0] < B[0] && X[1] >= A[1] && X[1] < B[1];
+    std::array<T, N> x;
+    for (int i = 0; i < N; i ++)
+      x[i] = X[i];
+    return contains(x);
   }
 
+  // bool contains(const T X[]) const {
+  //   return X[0] >= A[0] && X[0] < B[0] && X[1] >= A[1] && X[1] < B[1];
+  // }
+
   void update_centroid() {
-    C[0] = (A[0] + B[0]) / 2;
-    C[1] = (A[1] + B[1]) / 2;
+    for (int i = 0; i < N; i ++)
+      C[i] = (A[i] + B[i]) / 2;
+
+    // C[0] = (A[0] + B[0]) / 2;
+    // C[1] = (A[1] + B[1]) / 2;
   }
 
   void print() const {
-    fprintf(stderr, "A={%f, %f}, B={%f, %f}, centroid={%f, %f}\n", 
-        A[0], A[1], B[0], B[1], C[0], C[1]);
+    // fprintf(stderr, "A={%f, %f}, B={%f, %f}, centroid={%f, %f}\n", 
+    //     A[0], A[1], B[0], B[1], C[0], C[1]);
   }
 };
 
diff --git a/include/ftk/mesh/point_locator_2d_quad.hh b/include/ftk/mesh/point_locator_2d_quad.hh
index 67bef21f..d5a531bd 100644
--- a/include/ftk/mesh/point_locator_2d_quad.hh
+++ b/include/ftk/mesh/point_locator_2d_quad.hh
@@ -23,8 +23,8 @@ protected:
   struct quad_node {
     quad_node *parent = NULL;
     quad_node *children[4] = {NULL};
-    AABB<F> aabb;
-    std::vector<AABB<F>> elements; // valid only for leaf nodes
+    AABB<2, F> aabb;
+    std::vector<AABB<2, F>> elements; // valid only for leaf nodes
 
     ~quad_node();
     bool is_leaf() const { return elements.size() > 0; }
@@ -234,7 +234,7 @@ void point_locator_2d_quad<I, F>::initialize()
   root = new quad_node;
 
   // global bounds
-  AABB<F> &aabb = root->aabb;
+  AABB<2, F> &aabb = root->aabb;
   for (int i = 0; i < m2.n(0); i ++) {
     aabb.A[0] = std::min(aabb.A[0], coords[2*i]);
     aabb.A[1] = std::min(aabb.A[1], coords[2*i+1]);
@@ -244,7 +244,7 @@ void point_locator_2d_quad<I, F>::initialize()
   aabb.update_centroid();
   // aabb.print();
 
-  std::vector<AABB<F>> triangles(m2.n(2));
+  std::vector<AABB<2, F>> triangles(m2.n(2));
   for (int i=0; i<m2.n(2); i++) {
     const int i0 = conn[i*3], i1 = conn[i*3+1], i2 = conn[i*3+2];
     double x0 = coords[i0*2], x1 = coords[i1*2], x2 = coords[i2*2],