diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 652ca9f0f0db52fdb7dae21474fc1727ef26c745..2515ccdbf3379e971f505d516c757333b1f01d4d 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -105,6 +105,10 @@ namespace AMDiS {
     if (partStr == "zoltan")
       partitioner = new ZoltanPartitioner(&mpiComm);
 
+    tmp = 0;
+    GET_PARAMETER(0, name + "->box partitioning", "%d", &tmp);
+    partitioner->setBoxPartitioning(static_cast<bool>(tmp));
+
     TEST_EXIT(partitioner)("Could not create partitioner \"%s\"!\n", partStr.c_str());
   }
 
diff --git a/AMDiS/src/parallel/MeshPartitioner.cc b/AMDiS/src/parallel/MeshPartitioner.cc
index 5a469b1a928122a863d2367eec12cca106923384..fb832a2810764108a9b1d6538447e3c123237a7e 100644
--- a/AMDiS/src/parallel/MeshPartitioner.cc
+++ b/AMDiS/src/parallel/MeshPartitioner.cc
@@ -29,40 +29,99 @@ namespace AMDiS {
     // === Create initial partitioning of the AMDiS mesh. ===
 
     elementInRank.clear();
+    map<DofEdge, set<int> > vertexElements;
 
     TraverseStack stack;
     ElInfo *elInfo = 
       stack.traverseFirst(mesh, 0, 
-			  Mesh::CALL_EL_LEVEL | Mesh::FILL_NEIGH |  Mesh::FILL_BOUND);
+			  Mesh::CALL_EL_LEVEL | Mesh::FILL_NEIGH | Mesh::FILL_COORDS |  Mesh::FILL_BOUND);
     while (elInfo) {
-      int elIndex = elInfo->getElement()->getIndex();
-      int elInRank = std::min(elIndex / elPerRank, mpiSize - 1);
+      Element *el = elInfo->getElement();
+      int elIndex = el->getIndex();
 
-      if (elInRank == mpiRank)
-	elementInRank[elIndex] = true;
-      else
-	elementInRank[elIndex] = false;
-
-      partitionMap[elIndex] = elInRank;
-      
       for (int i = 0; i < mesh->getGeo(NEIGH); i++)
 	if (elInfo->getNeighbour(i) && elInfo->getBoundary(i) == INTERIOR) 
 	  elNeighbours[elIndex].push_back(elInfo->getNeighbour(i)->getIndex());
+
+      if (boxPartitioning) {
+	vertexElements[el->getEdge(0)].insert(elIndex);
+      } else {
+	int elInRank = std::min(elIndex / elPerRank, mpiSize - 1);
+	
+	elementInRank[elIndex] = (elInRank == mpiRank);
+	partitionMap[elIndex] = elInRank;	
+      }
+
       
       elInfo = stack.traverseNext(elInfo);
     }
+
+
+
+    if (boxPartitioning) {
+      TEST_EXIT(mesh->getDim() == 3)("Box partitioning only implemented for 3D!\n");
+
+      int boxCounter = 0;
+      for (map<DofEdge, set<int> >::iterator it = vertexElements.begin();
+	   it != vertexElements.end(); ++it) {
+	TEST_EXIT_DBG(it->second.size() == 6)("Should not happen!\n");
+	
+	boxSplitting[boxCounter] = it->second;
+
+	for (set<int>::iterator elIt = it->second.begin();
+	     elIt != it->second.end(); ++elIt)
+	  elInBox[*elIt] = boxCounter;
+
+	boxCounter++;
+      }
+
+      for (map<int, int>::iterator it = elInBox.begin(); it != elInBox.end(); ++it) {
+	int elBoxNo = it->second;
+
+	for (vector<int>::iterator neighIt = elNeighbours[it->first].begin();
+	     neighIt != elNeighbours[it->first].end(); ++neighIt) {
+	  int neighBoxNo = elInBox[*neighIt];
+
+	  if (elBoxNo != neighBoxNo)
+	    boxNeighbours[elBoxNo].insert(neighBoxNo);
+	}
+      }
+
+      MSG("Box partitioning with %d boxes enabled!\n", boxCounter);
+
+      int boxPerRank = boxCounter / mpiSize;
+
+      for (map<int, set<int> >::iterator it = boxSplitting.begin();
+	   it != boxSplitting.end(); ++it) {
+	int boxInRank = std::min(it->first / boxPerRank, mpiSize - 1);
+
+	for (set<int>::iterator elIt = it->second.begin();
+	     elIt != it->second.end(); ++elIt) {
+	  elementInRank[*elIt] = (boxInRank == mpiRank);
+	  partitionMap[*elIt] = boxInRank;
+	}
+      }
+    }
   }
 
 
   void MeshPartitioner::serialize(std::ostream &out)
   {
     SerUtil::serialize(out, elementInRank);
+    SerUtil::serialize(out, boxPartitioning);
+    SerUtil::serialize(out, boxSplitting);
+    SerUtil::serialize(out, boxNeighbours);
+    SerUtil::serialize(out, elInBox);
   }
 
 
   void MeshPartitioner::deserialize(std::istream &in)
   {
     SerUtil::deserialize(in, elementInRank);
+    SerUtil::deserialize(in, boxPartitioning);
+    SerUtil::deserialize(in, boxSplitting);
+    SerUtil::deserialize(in, boxNeighbours);
+    SerUtil::deserialize(in, elInBox);
   }
 
 }
diff --git a/AMDiS/src/parallel/MeshPartitioner.h b/AMDiS/src/parallel/MeshPartitioner.h
index 52701c7420c7f293520215cac03487d5519775f8..1bb8706e07a2360f187dc9a7e32e4b1c3212479e 100644
--- a/AMDiS/src/parallel/MeshPartitioner.h
+++ b/AMDiS/src/parallel/MeshPartitioner.h
@@ -24,6 +24,7 @@
 #define AMDIS_MESH_PARTITIONER_H
 
 #include <map>
+#include <set>
 #include "AMDiS_fwd.h"
 #include "Mesh.h"
 
@@ -47,6 +48,7 @@ namespace AMDiS {
     MeshPartitioner(MPI::Intracomm *comm)
       : mpiComm(comm),
 	mesh(NULL),
+	boxPartitioning(false),
 	mapLocalGlobal(NULL)
     {}
 
@@ -71,6 +73,11 @@ namespace AMDiS {
       mesh = m;
     }
 
+    void setBoxPartitioning(bool b)
+    {
+      boxPartitioning = b;
+    }
+
     void setLocalGlobalDofMap(map<DegreeOfFreedom, DegreeOfFreedom> *m)
     {
       mapLocalGlobal = m;
@@ -101,6 +108,14 @@ namespace AMDiS {
 
     Mesh *mesh;
 
+    bool boxPartitioning;
+
+    map<int, std::set<int> > boxSplitting;
+
+    map<int, std::set<int> > boxNeighbours;
+
+    map<int, int> elInBox;
+
     map<DegreeOfFreedom, DegreeOfFreedom> *mapLocalGlobal;
 
     map<int, vector<int> > elNeighbours;
diff --git a/AMDiS/src/parallel/ZoltanPartitioner.cc b/AMDiS/src/parallel/ZoltanPartitioner.cc
index c8c85fcb4ecbda12d0203ea9c338b4894236c262..e652e69aaf36daed3a1e30c2c0f62df198baa6c8 100644
--- a/AMDiS/src/parallel/ZoltanPartitioner.cc
+++ b/AMDiS/src/parallel/ZoltanPartitioner.cc
@@ -21,14 +21,7 @@ namespace AMDiS {
     : MeshPartitioner(comm),
       zoltan(*comm),
       elWeights(NULL)
-  {
-    zoltan.Set_Num_Obj_Fn(ZoltanFunctions::getNumObj, this);
-    zoltan.Set_Obj_List_Fn(ZoltanFunctions::getObjectList, this);
-    zoltan.Set_Num_Geom_Fn(ZoltanFunctions::getNumGeom, this);
-    zoltan.Set_Geom_Multi_Fn(ZoltanFunctions::getGeomMulti, this);
-    zoltan.Set_Num_Edges_Multi_Fn(ZoltanFunctions::getNumEdgeMulti, this);
-    zoltan.Set_Edge_List_Multi_Fn(ZoltanFunctions::getEdgeListMulti, this);
-  }
+  {}
 
 
   bool ZoltanPartitioner::partition(map<int, double> &weights,
@@ -36,6 +29,22 @@ namespace AMDiS {
   {
     FUNCNAME("ZoltanPartitioner::partition()");
 
+    if (!boxPartitioning) {
+      zoltan.Set_Num_Obj_Fn(ZoltanFunctions::getNumObj, this);
+      zoltan.Set_Obj_List_Fn(ZoltanFunctions::getObjectList, this);
+      zoltan.Set_Num_Geom_Fn(ZoltanFunctions::getNumGeom, this);
+      zoltan.Set_Geom_Multi_Fn(ZoltanFunctions::getGeomMulti, this);
+      zoltan.Set_Num_Edges_Multi_Fn(ZoltanFunctions::getNumEdgeMulti, this);
+      zoltan.Set_Edge_List_Multi_Fn(ZoltanFunctions::getEdgeListMulti, this);
+    } else {
+      zoltan.Set_Num_Obj_Fn(ZoltanFunctions::box_getNumObj, this);
+      zoltan.Set_Obj_List_Fn(ZoltanFunctions::box_getObjectList, this);
+      zoltan.Set_Num_Geom_Fn(ZoltanFunctions::box_getNumGeom, this);
+      zoltan.Set_Geom_Multi_Fn(ZoltanFunctions::box_getGeomMulti, this);
+      zoltan.Set_Num_Edges_Multi_Fn(ZoltanFunctions::box_getNumEdgeMulti, this);
+      zoltan.Set_Edge_List_Multi_Fn(ZoltanFunctions::box_getEdgeListMulti, this);
+    }
+
     if (mode != INITIAL)
       elWeights = &weights;
 
@@ -56,6 +65,8 @@ namespace AMDiS {
 
     if (mode == INITIAL) {
       zoltan.Set_Param("LB_APPROACH", "PARTITION");
+      if (boxPartitioning)
+	zoltan.Set_Param("LB_METHOD", "GRAPH");
     } else {
        zoltan.Set_Param("LB_APPROACH", "REPARTITION");
        zoltan.Set_Param("LB_METHOD", "GRAPH");
@@ -89,8 +100,17 @@ namespace AMDiS {
 	for (int i = 0; i < nImportEls; i++) {
 	  int recvElIndex = import_global_ids[i];
 	  
-	  elementInRank[recvElIndex] = true;
-	  recvElements[import_procs[i]].push_back(recvElIndex);
+	  if (!boxPartitioning) {
+	    elementInRank[recvElIndex] = true;
+	    recvElements[import_procs[i]].push_back(recvElIndex);
+	  } else {
+	    set<int>& box = boxSplitting[recvElIndex];
+
+	    for (set<int>::iterator it = box.begin(); it != box.end(); ++it) {
+	      elementInRank[*it] = true;
+	      recvElements[import_procs[i]].push_back(*it);
+	    }
+	  }
 	}
       }
       
@@ -98,8 +118,17 @@ namespace AMDiS {
 	for (int i = 0; i < nExportEls; i++) {
 	  int sendElIndex = export_global_ids[i];
 	  
-	  elementInRank[sendElIndex] = false;
-	  sendElements[export_procs[i]].push_back(sendElIndex);
+	  if (!boxPartitioning) {
+	    elementInRank[sendElIndex] = false;
+	    sendElements[export_procs[i]].push_back(sendElIndex);
+	  } else {
+	    set<int>& box = boxSplitting[sendElIndex];
+
+	    for (set<int>::iterator it = box.begin(); it != box.end(); ++it) {
+	      elementInRank[*it] = false;
+	      sendElements[export_procs[i]].push_back(*it);
+	    }
+	  }
 	}      
       }
 
@@ -167,8 +196,11 @@ namespace AMDiS {
 	partitionMap[partitionElements[elDist[i] + j]] = i;
   }
 
-  
 
+
+  // === Zoltan query function for "normal" element partitioning. ===
+
+  
   int ZoltanFunctions::getNumObj(void *data, int *ierr)
   {
     ZoltanPartitioner* zoltan = (ZoltanPartitioner*)data;
@@ -340,4 +372,166 @@ namespace AMDiS {
     *ierr = ZOLTAN_OK;
   }
 
+
+
+  // === Zoltan query functions for box partitioning. ===
+
+  int ZoltanFunctions::box_getNumObj(void *data, int *ierr)
+  {
+    FUNCNAME("ZoltanFunctions::box_getNumObj()");
+
+    ZoltanPartitioner* zoltan = (ZoltanPartitioner*)data;
+    map<int, bool>& elInRank = zoltan->getElementInRank();
+
+    int nObjects = 0;
+    for (map<int, bool>::iterator it = elInRank.begin(); it != elInRank.end(); ++it)
+      if (it->second)
+	nObjects++;
+
+    TEST_EXIT_DBG(nObjects % 6 == 0)("Should not happen!\n");
+
+    nObjects = nObjects / 6;
+
+    *ierr = ZOLTAN_OK;
+
+    return nObjects;
+  }
+
+
+  void ZoltanFunctions::box_getObjectList(void *data, 
+					  int sizeGid, 
+					  int sizeLid,
+					  ZOLTAN_ID_PTR globalId, 
+					  ZOLTAN_ID_PTR localId,
+					  int wgt_dim, 
+					  float *obj_wgts, 
+					  int *ierr)
+  {
+    ZoltanPartitioner* zoltan = (ZoltanPartitioner*)data;
+    map<int, bool>& elInRank = zoltan->getElementInRank();
+    set<int> boxes;
+    map<int, double> boxWeight;
+
+    for (map<int, bool>::iterator it = elInRank.begin(); it != elInRank.end(); ++it)
+      if (it->second) {
+	int boxId = zoltan->elInBox[it->first];
+	boxes.insert(boxId);
+	
+	if (zoltan->elWeights) {
+	  TEST_EXIT_DBG(zoltan->elWeights->count(it->first))
+	    ("No element weight for element index %d\n", it->first);
+	  
+	  boxWeight[boxId] += (*zoltan->elWeights)[it->first];
+	}
+      }
+
+    int localCounter = 0;
+    for (set<int>::iterator it = boxes.begin(); it != boxes.end(); ++it) {
+      globalId[localCounter] = *it;
+      localId[localCounter] = localCounter;
+	
+      if (wgt_dim > 0) {
+	TEST_EXIT_DBG(wgt_dim == 1)("Should not happen!\n");
+	
+	if (zoltan->elWeights)
+	  obj_wgts[localCounter] = static_cast<float>(boxWeight[*it]);
+	else
+	  obj_wgts[localCounter] = 1.0;	
+      }
+      
+      localCounter++;
+    }
+
+    *ierr = ZOLTAN_OK;
+  }
+
+
+  int ZoltanFunctions::box_getNumGeom(void *data, int *ierr)
+  {
+    ZoltanPartitioner* zoltan = (ZoltanPartitioner*)data;
+    *ierr = ZOLTAN_OK;
+
+    return zoltan->getMesh()->getGeo(WORLD);
+  }
+
+
+  void ZoltanFunctions::box_getGeomMulti(void *data, 
+					 int num_gid_entries, 
+					 int num_lid_entries, 
+					 int num_obj, 
+					 ZOLTAN_ID_PTR global_ids, 
+					 ZOLTAN_ID_PTR local_ids, 
+					 int num_dim, 
+					 double *geom_vec, 
+					 int *ierr)
+  {
+    FUNCNAME("ZoltanFunctions::box_getGeomMulti()");
+
+    ERROR_EXIT("Not yet supported!\n");
+  }
+
+
+  void ZoltanFunctions::box_getNumEdgeMulti(void *data, 
+					    int num_gid_entries, 
+					    int num_lid_entries, 
+					    int num_obj, 
+					    ZOLTAN_ID_PTR global_ids, 
+					    ZOLTAN_ID_PTR local_ids, 
+					    int *num_edges, 
+					    int *ierr)
+  {
+    FUNCNAME("ZoltanFunctions::box_getNumEdgeMulti()");
+
+    ZoltanPartitioner* zoltan = (ZoltanPartitioner*)data;
+
+    for (int i = 0; i < num_obj; i++) {
+      TEST_EXIT_DBG(zoltan->boxNeighbours.count(global_ids[i]))("Should not happen!\n");
+
+      num_edges[i] = zoltan->boxNeighbours[global_ids[i]].size();
+    }
+
+    *ierr = ZOLTAN_OK;
+  }
+
+
+  void ZoltanFunctions::box_getEdgeListMulti(void *data, 
+					     int num_gid_entries, 
+					     int num_lid_entries, 
+					     int num_obj, 
+					     ZOLTAN_ID_PTR global_ids, 
+					     ZOLTAN_ID_PTR local_ids, 
+					     int *num_edges, 
+					     ZOLTAN_ID_PTR nbor_global_id, 
+					     int *nbor_procs, 
+					     int wgt_dim, 
+					     float *ewgts, 
+					     int *ierr)
+  {
+    FUNCNAME("ZoltanFunctions::box_getEdgeListMulti()");
+
+    ZoltanPartitioner* zoltan = (ZoltanPartitioner*)data;
+
+    int c = 0;
+    for (int i = 0; i < num_obj; i++) {
+      TEST_EXIT_DBG(static_cast<int>(zoltan->boxNeighbours[global_ids[i]].size()) == 
+		    num_edges[i])
+	("Wrong sizes for global el index %d: %d %d\n",
+	 global_ids[i],
+	 zoltan->boxNeighbours[global_ids[i]].size(),
+	 num_edges[i]);
+
+      set<int> &neigh = zoltan->boxNeighbours[global_ids[i]];
+      for (set<int>::iterator it = neigh.begin(); it != neigh.end(); ++it) {
+	nbor_global_id[c] = *it;
+	nbor_procs[c] = zoltan->partitionMap[*(zoltan->boxSplitting[*it].begin())];
+	c++;
+      }
+    }
+
+    TEST_EXIT_DBG(wgt_dim == 0)("Edge weights are not handled yet!\n");
+
+    *ierr = ZOLTAN_OK;
+  }
+
+
 }
diff --git a/AMDiS/src/parallel/ZoltanPartitioner.h b/AMDiS/src/parallel/ZoltanPartitioner.h
index 909c2c8fdacab4e49eb47f4d83c1f4754111fab0..86262d6c86b67f3eb812d53824935b66311217d3 100644
--- a/AMDiS/src/parallel/ZoltanPartitioner.h
+++ b/AMDiS/src/parallel/ZoltanPartitioner.h
@@ -64,6 +64,8 @@ namespace AMDiS {
   class ZoltanFunctions
   {
   public:
+    // === Zoltan query function for "normal" element partitioning. ===
+
     static int getNumObj(void *data, int *ierr);   
 
     static void getObjectList(void *data, 
@@ -108,6 +110,55 @@ namespace AMDiS {
 				 int wgt_dim, 
 				 float *ewgts, 
 				 int *ierr);
+
+
+    // === Zoltan query functions for box partitioning. ===
+
+    static int box_getNumObj(void *data, int *ierr);   
+
+    static void box_getObjectList(void *data, 
+				  int sizeGid, 
+				  int sizeLid,
+				  ZOLTAN_ID_PTR globalId, 
+				  ZOLTAN_ID_PTR localId,
+				  int wgt_dim, 
+				  float *obj_wgts, 
+				  int *ierr);
+
+    static int box_getNumGeom(void *data, int *ierr);
+
+    static void box_getGeomMulti(void *data, 
+				 int num_gid_entries, 
+				 int num_lid_entries, 
+				 int num_obj, 
+				 ZOLTAN_ID_PTR global_ids, 
+				 ZOLTAN_ID_PTR local_ids, 
+				 int num_dim, 
+				 double *geom_vec, 
+				 int *ierr); 
+
+    static void box_getNumEdgeMulti(void *data, 
+				    int num_gid_entries, 
+				    int num_lid_entries, 
+				    int num_obj, 
+				    ZOLTAN_ID_PTR global_ids, 
+				    ZOLTAN_ID_PTR local_ids, 
+				    int *num_edges, 
+				    int *ierr);
+
+    static void box_getEdgeListMulti(void *data, 
+				     int num_gid_entries, 
+				     int num_lid_entries, 
+				     int num_obj, 
+				     ZOLTAN_ID_PTR global_ids, 
+				     ZOLTAN_ID_PTR local_ids, 
+				     int *num_edges, 
+				     ZOLTAN_ID_PTR nbor_global_id, 
+				     int *nbor_procs, 
+				     int wgt_dim, 
+				     float *ewgts, 
+				     int *ierr);
+
   };
 
 }