diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5a0257d..761b90c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.9.0] - 2022-10-09
+## Added
+- mesh2d as for diagnostic purposes
+- particle writer
+- docs
+- hypersonic test case
+- backup boundary module
+
+## [0.8.0] - 2022-10-03
+### Added
+- planar writer
+- planar diagnostic
+
 ## [0.7.0] - 2022-10-01
 ### Added
 - inflow boundary condition
diff --git a/doc/Boundary.ipynb b/doc/Boundary.ipynb
new file mode 100644
index 0000000..f2bc8b3
--- /dev/null
+++ b/doc/Boundary.ipynb
@@ -0,0 +1,167 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "d49edbe1-fe92-4a67-82e9-78a877139d5d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import sympy as sym"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7f598afc-282c-4bb5-8d47-80c77fcc1178",
+   "metadata": {},
+   "source": [
+    "**Boundary collision method**\n",
+    "\n",
+    "Given a plane $\\Pi$ and and a line $L$ as\n",
+    "\n",
+    "$$\n",
+    "\\Pi \\rightarrow (p - p_0) \\cdot \\vec{n}_p = 0 \\\\\n",
+    "L \\rightarrow l = l_0 + t \\cdot \\vec{n}_l\n",
+    "$$\n",
+    "\n",
+    "where $\\vec{n}_1$ is the normal vector of a plane given as\n",
+    "\n",
+    "$$\n",
+    "\\vec{n}_p = (p_1 - p_0) \\times (p_2 - p_0)\n",
+    "$$\n",
+    "\n",
+    "is a vector in the direction of the line as\n",
+    "\n",
+    "$$\n",
+    "\\vec{n}_l = l_1 - l_0.\n",
+    "$$\n",
+    "\n",
+    "The positons of the intersection point found by setting the point $p$ in the equation for the plane as being a point in the equation of the line:\n",
+    "\n",
+    "$$\n",
+    "(l - p_0) \\cdot \\vec{n}_p = 0 = (l_0 + t \\cdot \\vec{n}_l - p_0) \\cdot \\vec{n}_p \\implies t = - \\frac{(l_0 - p_0) \\cdot \\vec{n}_p}{\\vec{n}_p \\cdot \\vec{n}_l}\n",
+    "$$\n",
+    "\n",
+    "inserting the line equation the intersection point can be found as\n",
+    "\n",
+    "$$\n",
+    "l^* = l_0 + t \\cdot \\vec{n}_l.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "4abb0e41-5e15-4326-b663-df0b6b5cfe0b",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[0. 0. 0.]\n"
+     ]
+    }
+   ],
+   "source": [
+    "l1 = np.array([1.0, 1.0, 0.0])\n",
+    "l0 = np.array([-1.0, -1.0, 0.0])\n",
+    "\n",
+    "p0 = np.array([0.0, -1.0, -1.0])\n",
+    "p1 = np.array([0.0, -1.0, 1.0])\n",
+    "p2 = np.array([0.0, 1.0, 1.0])\n",
+    "\n",
+    "nl = l1 - l0\n",
+    "n_p = np.cross((p1 - p0), (p2 - p0))\n",
+    "\n",
+    "t = - ((l0 - p0).dot(n_p) / n_p.dot(nl))\n",
+    "\n",
+    "ls = l0 +t*nl\n",
+    "\n",
+    "print(ls)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d43f6e71-7324-4b21-9c0d-371fb76d01f6",
+   "metadata": {},
+   "source": [
+    "**Reflections**\n",
+    "\n",
+    "The reflected line can be found as\n",
+    "\n",
+    "$$\n",
+    "L_R \\rightarrow l_R = l^* + \\lambda \\cdot \\vec{n}_R\n",
+    "$$\n",
+    "\n",
+    "with \n",
+    "\n",
+    "$$\n",
+    "\\vec{n}_R = \\vec{n}_l - 2 \\frac{\\vec{n}_p \\cdot \\vec{n}_l}{|| \\vec{n}_p ||^2} \\cdot \\vec{n}_p.\n",
+    "$$\n",
+    "\n",
+    "The reflected point can now the found as \n",
+    "\n",
+    "$$\n",
+    "l_R = l^* + (1 - t) \\cdot \\vec{n}_R\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "91de06b7-7511-471d-9b64-f9c2a784d682",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([ 1., -1.,  0.])"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "nR = nl - 2*(n_p.dot(nl) / n_p.dot(n_p))*n_p\n",
+    "lR = ls - (1  -t )*nR\n",
+    "\n",
+    "lR"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a81e9e3a-d8d0-463f-a5e4-b8700e696875",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/dsmc/boundary.py b/dsmc/boundary.py
new file mode 100644
index 0000000..b2f979a
--- /dev/null
+++ b/dsmc/boundary.py
@@ -0,0 +1,180 @@
+import numpy as np
+from numba import njit
+from . import common as co
+
+@njit
+def _check_if_parallel(v1, v2, diff=1e-6):
+    n1 = np.linalg.norm(v1)
+    n2 = np.linalg.norm(v2)
+    
+    if n1 < 1.0e-13 or n2 < 1.0e-13:
+        return False
+    
+    V1 = np.copy(v1) / n1
+    V2 = np.copy(v2) / n2
+    
+    return V1.dot(V2) > diff
+
+@njit
+def _intersect(l0, l1, p0, p1, p2):
+    """
+    Args:
+        l0 : first position on line
+        l1 : second position on line
+        p0 : first position on plane
+        p1 : first position on plane
+        p2 : first position on plane
+        
+    Returns:
+        (intersected, n_l, n_r, t)
+    """
+    n_l = l1 - l0
+    n_p = np.cross((p1 - p0), (p2 - p1))
+    
+    if _check_if_parallel(n_l, n_p):
+        return (True, n_l, n_p, - ((l0 - p0).dot(n_p) / n_p.dot(n_l)))
+    else:
+        return (False, n_l, n_p, 0.0)
+
+@njit    
+def _calc_nr(n_l, n_p):
+    return n_l - 2.0 * (n_p.dot(n_l) / n_p.dot(n_p))*n_p
+
+@njit     
+def _reflect(vel, pos, pos_old, p0, p1, p2, domain):
+    intersected, n_l, n_p, t = _intersect(pos_old, pos, p0, p1, p2)
+    
+    if intersected and t < 1.0 and t > 0.0:
+        k = 1.0
+        p = pos_old + n_l*(t*k)
+        while not co.is_inside(p, domain):
+            k *= 0.9
+            p = pos_old + n_l*(t*k)
+        pos_old = p
+        n_r = _calc_nr(n_l, n_p)
+        pos = pos_old  + (1.0 - t)*n_r
+        vel = (np.linalg.norm(vel) / np.linalg.norm(n_r))*n_r
+
+    return (vel, pos, pos_old)
+
+@njit
+def _get_plane(domain, i, j):
+    if i == 0:
+        if j == 0:
+            p0 = np.array([domain[i][j], domain[1][0], domain[2][0]])
+            p1 = np.array([domain[i][j], domain[1][0], domain[2][1]])
+            p2 = np.array([domain[i][j], domain[1][1], domain[2][0]])
+        elif j == 1:
+            p0 = np.array([domain[i][j], domain[1][0], domain[2][0]])
+            p1 = np.array([domain[i][j], domain[1][1], domain[2][0]])
+            p2 = np.array([domain[i][j], domain[1][0], domain[2][1]])
+    elif i == 1:
+        if j == 0:
+            p0 = np.array([domain[0][0], domain[i][j], domain[2][0]])
+            p1 = np.array([domain[0][1], domain[i][j], domain[2][0]])
+            p2 = np.array([domain[0][0], domain[i][j], domain[2][1]])
+        if j == 1:
+            p0 = np.array([domain[0][0], domain[i][j], domain[2][0]])
+            p1 = np.array([domain[0][0], domain[i][j], domain[2][1]])
+            p2 = np.array([domain[0][1], domain[i][j], domain[2][0]])
+    elif i == 2:
+        if j == 0:
+            p0 = np.array([domain[0][0], domain[1][0], domain[i][j]])
+            p1 = np.array([domain[0][0], domain[1][1], domain[i][j]])
+            p2 = np.array([domain[0][1], domain[1][0], domain[i][j]])
+        if j == 1:
+            p0 = np.array([domain[0][0], domain[1][0], domain[i][j]])
+            p1 = np.array([domain[0][1], domain[1][0], domain[i][j]])
+            p2 = np.array([domain[0][0], domain[1][1], domain[i][j]])
+        
+    return (p0, p1, p2)
+
+@njit
+def _boundary(velocities, positions, old_positions, domain, boundary_conds):
+    kept_parts = np.ones(positions.shape[0], dtype=np.uint)
+    
+    for p in range(len(positions)):
+        counter = 0
+        while not co.is_inside(positions[p], domain) and kept_parts[p]:
+            if counter > 10:
+                kept_parts[p] = 0
+                break
+            
+            for i in range(3):
+                for j in range(2):
+                    p0, p1, p2 = _get_plane(domain, i, j)
+                    if boundary_conds[i][j] == 0:
+                        velocities[p], positions[p], old_positions[p] = _reflect(velocities[p], positions[p], old_positions[p], p0, p1, p2, domain)
+                        counter += 1
+                    elif boundary_conds[i][j] == 1 or boundary_conds[i][j] == 2:
+                        if _intersect(old_positions[p], positions[p], p0, p1, p2)[0]:
+                            kept_parts[p] = 0
+                        
+    N = int(sum(kept_parts))
+    p = 0
+    new_velocities = np.empty((N, 3))
+    new_positions = np.empty((N, 3))
+    new_old_positions = np.empty((N, 3))
+    
+    for i in range(positions.shape[0]):
+        if kept_parts[i] == 1:
+            new_velocities[p] = velocities[i]
+            new_positions[p] = positions[i]
+            new_old_positions[p] = old_positions[p]
+            p += 1
+        else:
+            continue
+
+    return (new_velocities, new_positions, new_old_positions)
+
+@njit
+def _get_boundary(boundary):
+    if boundary == "xmin":
+        return (0, 0)
+    elif boundary == "xmax":
+        return (0, 1)
+    elif boundary == "ymin":
+        return (1, 0)
+    elif boundary == "ymax":
+        return (1, 1)
+    elif boundary == "zmin":
+        return (2, 0)
+    elif boundary == "zmax":
+        return (2, 1)
+    
+@njit
+def _get_bc_type(bc_type):
+    if bc_type == "ela":
+        return 0
+    elif bc_type == "open":
+        return 1
+    elif bc_type == "inflow":
+        return 2
+    
+class Boundary:
+    def __init__(self):
+        self.T = np.ones((3, 2))*300.0
+        self.n = np.ones((3, 2))*1e+18
+        self.u = np.zeros((3, 2, 3))
+        self.boundary_conds = np.array([[0, 0], [0, 0], [0, 0]], dtype=np.uint) # 0 = ela, 1 = open, 2 = inflow 
+        self.domain = None
+        
+    def boundary(self, velocities, positions, old_positions):        
+        return _boundary(velocities, positions, old_positions, self.domain, self.boundary_conds)
+    
+    def set_bc_type(self, boundary, bc_type):
+        bound = _get_boundary(boundary)
+        bc = _get_bc_type(bc_type)
+        
+        self.boundary_conds[bound[0]][bound[1]] = bc
+        
+        print("boundary [" + boundary + "] set to [" + bc_type + "]")
+        
+    def set_bc_values(self, boundary, T, n, u):
+        i, j = _get_boundary(boundary)
+        
+        self.T[i][j] = T
+        self.n[i][j] = n
+        self.u[i][j] = u
+        
+        print("boundary [" + boundary + "] set to values T : {}, n : {}, u : {}".format(T, n, u))
\ No newline at end of file
diff --git a/dsmc/common.py b/dsmc/common.py
index 17bc81d..24e3fba 100644
--- a/dsmc/common.py
+++ b/dsmc/common.py
@@ -1,4 +1,5 @@
 from numba import cfunc, njit
+import numpy.typing as npt
 
 @cfunc("double(double[::3, ::2])")
 def get_V(a):
@@ -8,4 +9,12 @@ def get_V(a):
 def swap(arr, pos1, pos2):
     temp = arr[pos2]
     arr[pos2] = arr[pos1]
-    arr[pos1] = temp
\ No newline at end of file
+    arr[pos1] = temp
+    
+@njit
+def is_inside(position : npt.NDArray, box : npt.NDArray) -> bool:
+    a : bool = position[0] >= box[0][0] and position[0] <= box[0][1]
+    b : bool = position[1] >= box[1][0] and position[1] <= box[1][1]
+    c : bool = position[2] >= box[2][0] and position[2] <= box[2][1]
+
+    return a and b and c
\ No newline at end of file
diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py
index 97a7e58..6f5a8ac 100644
--- a/dsmc/diagnostics.py
+++ b/dsmc/diagnostics.py
@@ -1,64 +1,33 @@
 import numpy as np
-from numba import njit
 from . import particles as prt
 from . import octree as oc
-from . import common as com
-
-@njit
-def _in_2d_box(pos, box):
-    a = pos[0] >= box[0][0] and pos[0] <= box[0][1]
-    b = pos[1] >= box[1][0] and pos[1] <= box[1][1]
-    return a and b
-
-@njit
-def _set_up_2d_boxes(x0, x1, y0, y1, Nx, Ny):
-    dx = (x1 - x0) / Nx
-    dy = (y1 - y0) / Ny
-    boxes = np.empty((Nx*Ny, 2, 2))
-    y = y0
-    
-    for i in range(Ny):
-        x = x0
-        for j in range(Nx):
-            k = j + i*Nx
-            boxes[k][0][0] = x
-            boxes[k][0][1] = x + dx
-            boxes[k][1][0] = y
-            boxes[k][1][1] = y + dy
-            x += dx
-        y += dy
-    
-    return boxes
-
-@njit
-def _sort_2d(positions, x0, x1, y0, y1, Nx, Ny):
-    boxes = _set_up_2d_boxes(x0, x1, y0, y1, Nx, Ny) # [] : bix id, [][] : axis, [][][] : 0 min , 1 max
-    permutations = [i for i in range(len(positions))]
-    leafs = np.zeros((Nx * Ny, 2), dtype=np.int_) # leaf[0] : ofseat, leaf[1] : n-parts
-    
-    for i in range(len(leafs)):
-        leafs[i][0] = leafs[i - 1][0] + leafs[i - 1][1] if i > 0 else 0
-        leafs[i][1] = 0
-        runner = leafs[i][0]
-        for j in range(leafs[i][0], len(positions)):
-            p = permutations[j]
-            if _in_2d_box(positions[p], boxes[i]):
-                com.swap(permutations, j, runner)
-                runner += 1
-                leafs[i][1] += 1
-
-    return (permutations, leafs, boxes)
+#from . import common as com
 
 class Values:
     def __init__(self):
         self.number_elements = 0
+        
+    def calc_vals(self, positions, velocities, ids, box):
+        #V = com.get_V(box)
+        self.number_elements = len(ids)
 
-def analyse_2d(positions, x0, x1, y0, y1, Nx, Ny):
-    permutations, leafs, boxes = _sort_2d(positions, x0, x1, y0, y1, Nx, Ny)
-    values = [Values() for _ in range(len(leafs))]
+
+def analyse_2d(positions, velocities, mesh2d, h):
+    boxes = np.empty((len(mesh2d.cells), 3, 2))
+    values = [Values() for _ in range(len(mesh2d.cells))]
     
-    for i in range(len(leafs)):
-        values[i].number_elements = leafs[i][1]
+    for y in range(mesh2d.n_cells2):
+        for x in range(mesh2d.n_cells1):
+            N = x + y * mesh2d.n_cells1
+            
+            boxes[N][0][0] = mesh2d.min1 + mesh2d.cell_size1*x
+            boxes[N][0][1] = mesh2d.min1 + mesh2d.cell_size1*(x + 1)
+            boxes[N][1][0] = mesh2d.min2 + mesh2d.cell_size2*y
+            boxes[N][1][1] = mesh2d.min2 + mesh2d.cell_size2*(y + 1)
+            boxes[N][2][0] = 0.0
+            boxes[N][2][1] = h
+            
+            values[N].calc_vals(positions, velocities, mesh2d.cells[N], boxes[N])
         
     return (boxes, values)
                     
diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py
index 4b7d212..01209ae 100644
--- a/dsmc/dsmc.py
+++ b/dsmc/dsmc.py
@@ -7,12 +7,13 @@
 
 @njit
 def _push(velocities, positions, dt):
+    old_positions = np.copy(positions)
     for p in prange(len(positions)):
         positions[p] = positions[p] + velocities[p]*dt
-    return positions
+    return (velocities, positions, old_positions)
 
 @njit
-def _boundary(velocities, positions, domain, boundary_conds):
+def _boundary(velocities, positions, old_positions, domain, boundary_conds):
     kept_parts = np.ones(positions.shape[0], dtype=np.uint)
     
     for p in prange(len(positions)):
@@ -20,31 +21,35 @@ def _boundary(velocities, positions, domain, boundary_conds):
             for i in range(3):
                 if positions[p][i] < domain[i][0]:
                     if boundary_conds[i][0] == 0:
+                        old_positions[p][i] = positions[p][i]
                         positions[p][i] = 2.0 * domain[i][0] - positions[p][i]
                         velocities[p][i] *= -1.0
                     elif boundary_conds[i][0] == 1 or boundary_conds[i][0] == 2:
                         kept_parts[p] = 0
                 if positions[p][i] > domain[i][1]:
                     if boundary_conds[i][1] == 0:
+                        old_positions[p][i] = positions[p][i]
                         positions[p][i] = 2.0 * domain[i][1] - positions[p][i]
                         velocities[p][i] *= -1.0
                     elif boundary_conds[i][1] == 1 or boundary_conds[i][0] == 2:
                         kept_parts[p] = 0
                         
-    N = sum(kept_parts)
+    N = int(sum(kept_parts))
     p = 0
     new_velocities = np.empty((N, 3))
     new_positions = np.empty((N, 3))
+    new_old_positions = np.empty((N, 3))
     
     for i in range(positions.shape[0]):
         if kept_parts[i] == 1:
             new_velocities[p] = velocities[i]
             new_positions[p] = positions[i]
+            new_old_positions[p] = old_positions[p]
             p += 1
         else:
             continue
 
-    return (new_velocities, new_positions)
+    return (new_velocities, new_positions, new_old_positions)
 
 @njit
 def _check_positions(velocities, positions, old_positions, domain):
@@ -54,6 +59,31 @@ def _check_positions(velocities, positions, old_positions, domain):
         if (not oc._is_inside(positions[i], domain)) and (not oc._is_inside(old_positions[i], domain)):
             kept_parts[i] = 0
     
+    N = sum(kept_parts)
+    p = 0
+    new_velocities = np.empty((N, 3))
+    new_positions = np.empty((N, 3))
+    new_old_positions = np.empty((N, 3))
+    
+    for i in prange(positions.shape[0]):
+        if kept_parts[i] == 1:
+            new_velocities[p] = velocities[i]
+            new_positions[p] = positions[i]
+            new_old_positions[p] = old_positions[i]
+            p += 1
+        else:
+            continue
+
+    return (new_velocities, new_positions, new_old_positions)
+
+@njit
+def _check_created_particles(velocities, positions, obj):
+    kept_parts = np.ones(positions.shape[0], dtype=np.uint)
+    
+    for i in prange(positions.shape[0]):
+        if oc._is_inside(positions[i], obj):
+            kept_parts[i] = 0
+    
     N = sum(kept_parts)
     p = 0
     new_velocities = np.empty((N, 3))
@@ -69,6 +99,23 @@ def _check_positions(velocities, positions, old_positions, domain):
 
     return (new_velocities, new_positions)
 
+@njit
+def _object(velocities, positions, old_positions, coll_obj):
+    for p in range(positions.shape[0]):
+        if oc._is_inside(positions[p], coll_obj):
+            for i in range(3):
+                if (old_positions[p][i] < coll_obj[i][0]):
+                    old_positions[p][i] = positions[p][i]
+                    positions[p][i] = 2.0 * coll_obj[i][0] - positions[p][i]
+                    velocities[p][i] *= -1.0
+                
+                if (old_positions[p][i] > coll_obj[i][1]):
+                    old_positions[p][i] = positions[p][i]
+                    positions[p][i] = 2.0 * coll_obj[i][1] - positions[p][i]
+                    velocities[p][i] *= -1.0
+                    
+    return (velocities, positions, old_positions)
+
 @njit
 def _calc_prob(rel_vel : float, sigma_T : float, Vc : float, dt : float, w : float, N : int) -> np.single:
     """
@@ -179,12 +226,13 @@ def clear(self):
         self.boundary = Boundary()
         self.sigma_T = 3.631681e-19
         self.mass = None
+        self.objects = []
 
     def advance(self, dt, collisions=True, octree=True):
         if self.domain is None:
             raise Exception("simulation domain not defined")
         if self.particles.N == 0:
-            raise Exception("no particles in domain")
+            print("warning: no particles in domain")
         if self.w == None:
             raise Exception("particle weight not set")
             
@@ -193,29 +241,39 @@ def advance(self, dt, collisions=True, octree=True):
                 if self.boundary_conds[i][j] == 2:
                     self.particles.inflow(self.mass, self.boundary.T[i][j], self.boundary.u[i][j], self.boundary.n[i][j], self.w, dt, self.domain, i, j)
 
+        velocities, positions, old_positions = _push(self.particles.Vel, self.particles.Pos, dt)
+        velocities, positions, old_positions = _check_positions(velocities, positions, old_positions, self.domain)
+        
+        
+        for obj in self.objects:
+            velocities, positions, old_positions  = _object(velocities, positions, old_positions, obj)
+            
+        velocities, positions, old_positions = _boundary(velocities, positions, old_positions, self.domain, self.boundary_conds)
+            
         if octree:
-            self.octree.build(self.particles.Pos)
+            self.octree.build(positions)
         if collisions and octree:
-            self.particles.VelPos = (self._update_velocities(dt), self.particles.Pos)
-        old_positions = np.copy(self.particles.Pos)
-        positions = _push(self.particles.Vel, self.particles.Pos, dt)
-        self.particles.VelPos = _check_positions(self.particles.Vel, positions, old_positions, self.domain)
-        self.particles.VelPos = _boundary(self.particles.Vel, self.particles.Pos, self.domain, self.boundary_conds)
+            velocities = self._update_velocities(dt, velocities)
+             
+        self.particles.VelPos = (velocities, positions)
 
-    def _update_velocities(self, dt):
+    def _update_velocities(self, dt, velocities):
         Nleafs : int = len(self.octree.leafs)
         elem_offsets : np.ndarray = np.array([leaf.elem_offset for leaf in self.octree.leafs], dtype=int)
         number_elements : np.ndarray = np.array([leaf.number_elements for leaf in self.octree.leafs], dtype=int)
         number_children : np.ndarray = np.array([leaf.number_children for leaf in self.octree.leafs], dtype=int)
         cell_boxes : np.ndarray = np.array([box for box in self.octree.cell_boxes])
 
-        return _update_vels(self.octree.permutations, self.particles.Vel, self.mass, self.sigma_T, dt, self.w, elem_offsets, number_elements, number_children, cell_boxes, Nleafs)
+        return _update_vels(self.octree.permutations, velocities, self.mass, self.sigma_T, dt, self.w, elem_offsets, number_elements, number_children, cell_boxes, Nleafs)
 
     def create_particles(self, box, T, n, u = np.zeros(3)):
         box = np.array(box)
         N = int(round(com.get_V(box) * n / self.w))
         print("creating {} particles".format(N))
         self.particles.create_particles(box, self.mass, T, N, u)
+        
+        for obj in self.objects:
+            self.particles.VelPos = _check_created_particles(self.particles.Vel, self.particles.Pos, obj)
 
         print("now containing {} particles, {} total".format(N, self.particles.N))
 
@@ -245,3 +303,6 @@ def set_bc_values(self, boundary, T, n, u):
         self.boundary.u[i][j] = u
         
         print("boundary [" + boundary + "] set to values T : {}, n : {}, u : {}".format(T, n, u))
+
+    def add_object(self, coll_object):
+        self.objects.append(np.array(coll_object))
\ No newline at end of file
diff --git a/dsmc/mesh/__init__.py b/dsmc/mesh/__init__.py
new file mode 100644
index 0000000..633f866
--- /dev/null
+++ b/dsmc/mesh/__init__.py
@@ -0,0 +1,2 @@
+# -*- coding: utf-8 -*-
+
diff --git a/dsmc/mesh/mesh2d.py b/dsmc/mesh/mesh2d.py
new file mode 100644
index 0000000..a409084
--- /dev/null
+++ b/dsmc/mesh/mesh2d.py
@@ -0,0 +1,82 @@
+from enum import Enum
+import math
+from numba import njit
+import numpy as np
+
+@njit
+def _get_cell_id(val1, val2, n_cells1, n_cells2, min1, min2, cell_size1, cell_size2):
+    if (val1 < min1):
+        return (False, 0)
+    elif (val1 > (min1 + n_cells1 * cell_size1)):
+        return (False, 0)
+    else:
+        cell_id1 = math.floor((val1 - min1) / cell_size1)
+
+    if (val2 < min2):
+        return (False, 0)
+    elif (val2 > (min2 + n_cells2 * cell_size2)):
+        return (False, 0)
+    else:
+        cell_id2 = math.floor((val2 - min2) / cell_size2)
+
+    return (True, cell_id2 * n_cells1 + cell_id1)
+
+@njit
+def _sort(values1, values2, n_cells1, n_cells2, min1, min2, cell_size1, cell_size2):
+    inside = np.empty((len(values1)), np.bool_)
+    ids = np.empty((len(values1)), np.int_)
+    
+    for i in range(len(values1)):
+        inside[i], ids[i] = _get_cell_id(values1[i], values2[i], n_cells1, n_cells2, min1, min2, cell_size1, cell_size2)
+        
+    return (inside, ids)
+
+
+class Plane(Enum):
+    XY = 0
+    YZ = 1
+    XZ = 2
+
+class Mesh2d:
+    def __init__(self):
+        self.clear()
+        
+    def clear(self):
+        self.cell_size1 = 1.0
+        self.cell_size2 = 1.0
+        self.min1 = 0.0
+        self.min2 = 0.0
+        self.n_cells1 = 1
+        self.n_cells2 = 1
+        self.plane = Plane.XY
+        self.cells = None # this holds an 2d array, cells[] : cells, cells[][] position ids in cell
+        
+    def sort(self, position):
+        self.cells = [[] for _ in range(self.n_cells1 * self.n_cells2)]
+        
+        match self.plane:
+            case Plane.XY:
+                positions1 = np.array([position[i][0] for i in range(len(position))])
+                positions2 = np.array([position[i][1] for i in range(len(position))])
+            case Plane.YZ:
+                positions1 = np.array([position[i][1] for i in range(len(position))])
+                positions2 = np.array([position[i][2] for i in range(len(position))])
+            case Plane.XZ:
+                positions1 = np.array([position[i][0] for i in range(len(position))])
+                positions2 = np.array([position[i][2] for i in range(len(position))])
+            
+        inside, cell_ids =  _sort(positions1, positions2, self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size1, self.cell_size2)
+        sorted_ids = np.argsort(cell_ids)
+        
+        for idx in sorted_ids:
+            if inside[idx]:
+                self.cells[cell_ids[idx]].append(idx)
+        
+    def get_cell_id(self, position):
+        match self.plane:
+            case Plane.XY:
+                return _get_cell_id(position[0], position[1], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size1, self.cell_size2)
+            case Plane.YZ:
+                return _get_cell_id(position[1], position[2], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size1, self.cell_size2)
+            case Plane.XZ:
+                return _get_cell_id(position[0], position[2], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size1, self.cell_size2)
\ No newline at end of file
diff --git a/dsmc/particles.py b/dsmc/particles.py
index 421fbff..4516065 100644
--- a/dsmc/particles.py
+++ b/dsmc/particles.py
@@ -83,7 +83,7 @@ def calc_temperature(velocities, mass):
     return tot_e / ((3.0/2.0) * len(velocities) * kb)
 
 @njit
-def calc_positions(x, y, z, N):
+def calc_positions(X, N):
     """
     Parameters
     ----------
@@ -97,9 +97,8 @@ def calc_positions(x, y, z, N):
     positions = np.empty((N, 3))
 
     for i in range(N):
-        positions[i][0] = x[0] + np.random.random() * (x[1] - x[0])
-        positions[i][1] = y[0] + np.random.random() * (y[1] - y[0])
-        positions[i][2] = z[0] + np.random.random() * (z[1] - z[0])
+        for j in range(3):
+            positions[i][j] = X[j][0] + np.random.random() * (X[j][1] - X[j][0])
 
     return positions
 
@@ -136,11 +135,11 @@ def VelPos(self, vel_pos):
     def create_particles(self, X, mass, T, N, u = np.zeros(3)):    
         if self._N == 0:
             self._velocities = get_velocities(T, mass, N, u)
-            self._positions = calc_positions(X[0], X[1], X[2], N)
+            self._positions = calc_positions(X, N)
             self._N  = N
         else:
             self._velocities = np.concatenate((self._velocities, get_velocities(T, mass, N, u)))
-            self._positions = np.concatenate((self._positions, calc_positions(X[0], X[1], X[2], N)))
+            self._positions = np.concatenate((self._positions, calc_positions(X, N)))
             self._N  += N
             
 
diff --git a/dsmc/writer/__init__.py b/dsmc/writer/__init__.py
index 6847e5f..a53d8ad 100644
--- a/dsmc/writer/__init__.py
+++ b/dsmc/writer/__init__.py
@@ -1,2 +1,3 @@
 from .octree import write_buttom_leafs
 from .planar import write as write_planar
+from .particles import write_positions
diff --git a/dsmc/writer/particles.py b/dsmc/writer/particles.py
new file mode 100644
index 0000000..b6f2c5d
--- /dev/null
+++ b/dsmc/writer/particles.py
@@ -0,0 +1,5 @@
+def write_positions(particles, file_name="positions.csv"):
+    with open(file_name, "w") as file:
+        file.write("x, y, z\n")
+        for pos in particles.Pos:
+            file.write("{}, {}, {}\n".format(pos[0], pos[1], pos[2]))
\ No newline at end of file
diff --git a/examples/hypersonic_flow.py b/examples/hypersonic_flow.py
new file mode 100644
index 0000000..d084f39
--- /dev/null
+++ b/examples/hypersonic_flow.py
@@ -0,0 +1,90 @@
+import dsmc
+import dsmc.writer as wrt
+import dsmc.diagnostics as dia
+import dsmc.mesh.mesh2d as msh2d
+import dsmc.octree as oc
+import time
+import numpy as np
+
+if __name__ == '__main__':
+    # general parameters
+    solver = dsmc.DSMC()
+    domain = [(-3.0, 3.0), (-1.5, 1.5), (-0.025, 0.025)]
+    obj = [(-0.25, 0.25), (-0.25, 0.25), (-0.5, 0.5)]
+    dt = 1e-6
+    w = 0.2 * 2.3e+15
+    mass = 6.6422e-26
+    T =  273.0
+    n = 2.6e+19
+    u = np.array([0.0, -3043.0, 0.0])
+    niter = 1000
+    niter_sample = 100
+    
+    # set up mesh2
+    mesh = msh2d.Mesh2d()
+    
+    mesh.n_cells1 = 100
+    mesh.n_cells2 = 50
+    mesh.min1 = domain[0][0]
+    mesh.min2 = domain[1][0]
+    mesh.cell_size1 = 0.06
+    mesh.cell_size2 = 0.06
+    
+    h = domain[2][1] - domain[2][0]
+    
+    # setup solver
+    solver.set_domain(domain)
+    solver.set_weight(w)
+    solver.set_mass(mass)
+    
+    # set up boundary
+    solver.set_bc_type("xmin", "inflow")
+    solver.set_bc_type("xmax", "inflow")
+    solver.set_bc_type("ymax", "inflow")
+    
+    solver.set_bc_type("ymin", "open")
+    
+    solver.set_bc_values("xmin", T, n, u)
+    solver.set_bc_values("xmax", T, n, u)
+    solver.set_bc_values("ymax", T, n, u)
+    
+    # set object
+    solver.add_object(obj)
+    
+    # create particles
+    solver.create_particles(domain, T, n, u)
+    
+    # trees
+    tree_inner = oc.Octree()
+    tree_outer = oc.Octree()
+ 
+    inner_pos = np.array([[-0.25, -0.25, -0.025], [0.25, 0.25, 0.025]])
+    outer_pos = np.array([[-3.0, -1.5, -0.025], [3.0, 1.5, 0.025]])
+    
+    tree_inner.build(inner_pos)
+    tree_outer.build(outer_pos)
+  
+    wrt.write_buttom_leafs(tree_inner, "innter.vtu")
+    wrt.write_buttom_leafs(tree_outer, "outer.vtu")
+    
+    # start timing
+    start_time = time.time()
+    
+    for it in range(niter):
+        print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True)
+        solver.advance(dt)
+            
+    for it in range(niter_sample):
+        print("iteration {:4}/{}".format(it + 1, niter_sample), end="\r", flush=True)
+        solver.advance(dt)
+        
+        if (it + 1)%10 == 0:
+            mesh.sort(solver.particles.Pos)
+            boxes, values = dia.analyse_2d(solver.particles.Pos, solver.particles.Vel, mesh, h)
+            wrt.write_planar(boxes, values, "hypersonic_{}.vtu".format(it))
+        
+    wrt.write_buttom_leafs(solver.octree)
+        
+    print("")
+    print("--- %s seconds ---" % (time.time() - start_time))
+    print('done')
\ No newline at end of file
diff --git a/examples/hypersonic_flow_mini.py b/examples/hypersonic_flow_mini.py
new file mode 100644
index 0000000..109eb1c
--- /dev/null
+++ b/examples/hypersonic_flow_mini.py
@@ -0,0 +1,70 @@
+import dsmc
+import dsmc.writer as wrt
+import dsmc.octree as oc
+import time
+import numpy as np
+
+if __name__ == '__main__':
+    # general parameters
+    solver = dsmc.DSMC()
+    domain = [(-3.0, 3.0), (-1.5, 1.5), (-0.025, 0.025)]
+    obj = [(-0.25, 0.25), (-0.25, 0.25), (-0.5, 0.5)]
+    dt = 1e-6
+    w = 0.2 * 2.3e+15
+    mass = 6.6422e-26
+    T =  273.0
+    n = 2.6e+19
+    u = np.array([0.0, -3043.0, 0.0])
+    niter = 500
+    
+    h = domain[2][1] - domain[2][0]
+    
+    # setup solver
+    solver.set_domain(domain)
+    solver.set_weight(w)
+    solver.set_mass(mass)
+    
+    # set object
+    solver.add_object(obj)
+    
+    # trees
+    tree_inner = oc.Octree()
+    tree_outer = oc.Octree()
+    
+    inner_pos = np.array([[-0.25, -0.25, -0.025], [0.25, 0.25, 0.025]])
+    outer_pos = np.array([[-3.0, -1.5, -0.025], [3.0, 1.5, 0.025]])
+    
+    tree_inner.build(inner_pos)
+    tree_outer.build(outer_pos)
+    
+    # create particles
+    positions = np.zeros((2, 3))
+    velocities = np.zeros((2, 3))
+    
+    positions[0][0] = -0.01
+    positions[0][1] = 1.4
+    positions[0][2] = -0.01
+    
+    positions[1][0] = 0.01
+    positions[1][1] = 1.3
+    positions[1][2] = 0.01
+    
+    velocities[0][1] = -3043.0
+    velocities[1][1] = -3043.0
+    
+    solver.particles.VelPos = (velocities, positions)
+    
+    # start timing
+    start_time = time.time()
+    
+    for it in range(niter):
+        print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True)
+        solver.advance(dt)
+        wrt.write_positions(solver.particles, "pos_{}.csv".format(it))
+        
+    wrt.write_buttom_leafs(tree_inner, "innter.vtu")
+    wrt.write_buttom_leafs(tree_outer, "outer.vtu")
+        
+    print("")
+    print("--- %s seconds ---" % (time.time() - start_time))
+    print('done')
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 0b8444c..800507d 100644
--- a/setup.py
+++ b/setup.py
@@ -1,10 +1,10 @@
 from setuptools import setup
 setup(
     name='dsmc',
-    version='0.7.0',
+    version='0.9.0',
     author='Leo Basov',
     python_requires='>=3.10, <4',
-    packages=["dsmc", "dsmc.writer"],
+    packages=["dsmc", "dsmc.writer", "dsmc.mesh"],
     install_requires=[
         'numpy',
         'llvmlite',
diff --git a/tests/domain/boundary.py b/tests/domain/boundary.py
new file mode 100644
index 0000000..5d70ba9
--- /dev/null
+++ b/tests/domain/boundary.py
@@ -0,0 +1,93 @@
+import dsmc.dsmc as ds
+import dsmc.particles as prt
+import dsmc.octree as oc
+import dsmc.writer as wrt
+import dsmc.boundary as bo
+import dsmc.common as co
+import time
+import numpy as np
+
+def create_pos_and_vels():
+    positions = np.zeros((6, 3))
+    velocitiies = np.zeros((6, 3))
+    
+    # x
+    positions[0][0] = -0.75
+    velocitiies[0][0] = 100.0
+    
+    positions[1][0] = 0.75
+    velocitiies[1][0] = -100.0
+    
+    # y
+    positions[2][1] = -0.75
+    velocitiies[2][1] = 100.0
+    
+    positions[3][1] = 0.75
+    velocitiies[3][1] = -100.0
+    
+    # z
+    positions[4][2] = -0.75
+    velocitiies[4][2] = 100.0
+    
+    positions[5][2] = 0.75
+    velocitiies[5][2] = -100.0
+    
+    return (velocitiies, positions)
+
+if __name__ == '__main__':
+    # general parameters
+    boundary = bo.Boundary()
+    particles = prt.Particles()
+    domain = np.array([(-1.0, 1.0), (-1.0, 1.0), (-1.0, 1.0)])
+    dt = 0.0001
+    w = 2.3e+16
+    mass = 6.6422e-26
+    T =  273.0
+    n = 2.6e+19
+    N = int(co.get_V(domain)*n/w)
+    niter = 1000
+    
+    # trees
+    tree_outer = oc.Octree()
+    outer_pos = np.array([[-1.0, -1.0, -1.0], [1.0, 1.0, 1.0]])
+    tree_outer.build(outer_pos)
+    #wrt.write_buttom_leafs(tree_outer, "outer_box.vtu")
+    
+    # setup solver
+    boundary.domain = domain
+    
+    # create particles
+    particles.create_particles(domain, mass, T, N)
+    
+    velocities, positions = particles.VelPos
+        
+    #particles.VelPos = create_pos_and_vels()
+    
+    # start timing
+    start_time = time.time()
+    
+    wrt.write_positions(particles, "pos_{}.csv".format(0))
+    
+    E0 = 0.0
+    
+    for vel in particles.Vel:
+        E0 += vel.dot(vel)
+    
+    for it in range(niter):
+        E = 0.0
+        
+        velocities, positions, old_positions = ds._push(particles.Vel, particles.Pos, dt)
+        velocities, positions, old_positions = boundary.boundary(velocities, positions, old_positions)
+        
+        particles.VelPos = (velocities, positions)
+        wrt.write_positions(particles, "pos_{}.csv".format(it + 1))
+        
+        for vel in particles.Vel:
+            E += vel.dot(vel)
+        
+        print("iteration {:4}/{}, N particles {}/{}, Efrac {}".format(it + 1, niter, particles.N, N, E/E0), end="\r", flush=True)
+    
+
+    print("")
+    print("--- %s seconds ---" % (time.time() - start_time))
+    print('done')
\ No newline at end of file
diff --git a/tests/domain/domain_obj.py b/tests/domain/domain_obj.py
new file mode 100644
index 0000000..22c8ad1
--- /dev/null
+++ b/tests/domain/domain_obj.py
@@ -0,0 +1,86 @@
+import dsmc
+import dsmc.octree as oc
+import dsmc.writer as wrt
+import time
+import numpy as np
+
+def create_pos_and_vels():
+    positions = np.zeros((6, 3))
+    velocitiies = np.zeros((6, 3))
+    
+    # x
+    positions[0][0] = -0.75
+    velocitiies[0][0] = 1.0
+    
+    positions[1][0] = 0.75
+    velocitiies[1][0] = -1.0
+    
+    # y
+    positions[2][1] = -0.75
+    velocitiies[2][1] = 1.0
+    
+    positions[3][1] = 0.75
+    velocitiies[3][1] = -1.0
+    
+    # z
+    positions[4][2] = -0.75
+    velocitiies[4][2] = 1.0
+    
+    positions[5][2] = 0.75
+    velocitiies[5][2] = -1.0
+    
+    return (velocitiies, positions)
+
+def write_partices(positions, it):
+    with open("particles_{}.csv".format(it), "w") as file:
+        file.write("x, y, z\n")
+        for pos in positions:
+            file.write("{}, {}, {}\n".format(pos[0], pos[1], pos[2]))
+
+if __name__ == '__main__':
+    # general parameters
+    solver = dsmc.DSMC()
+    domain = [(-1.0, 1.0), (-1.0, 1.0), (-1.0, 1.0)]
+    obj = [(-0.5, 0.5), (-0.5, 0.5), (-0.5, 0.5)]
+    dt = 0.01
+    w = 2.3e+14
+    mass = 6.6422e-26
+    T =  273.0
+    n = 2.6e+19
+    niter = 200
+    
+    # trees
+    tree_inner = oc.Octree()
+    tree_outer = oc.Octree()
+    
+    inner_pos = np.array([[-0.5, -0.5, -0.5], [0.5, 0.5, 0.5]])
+    outer_pos = np.array([[-1.0, -1.0, -1.0], [1.0, 1.0, 1.0]])
+    
+    tree_inner.build(inner_pos)
+    tree_outer.build(outer_pos)
+    
+    # setup solver
+    solver.set_domain(domain)
+    solver.set_weight(w)
+    solver.set_mass(mass)
+    
+    # set object
+    solver.add_object(obj)
+    
+    # create particles
+    solver.particles.VelPos = create_pos_and_vels()
+    
+    # start timing
+    start_time = time.time()
+    
+    for it in range(niter):
+        print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True)
+        solver.advance(dt, octree=False)
+        write_partices(solver.particles.Pos, it)
+        
+    wrt.write_buttom_leafs(tree_inner, "inner_box.vtu")
+    wrt.write_buttom_leafs(tree_outer, "outer_box.vtu")
+
+    print("")
+    print("--- %s seconds ---" % (time.time() - start_time))
+    print('done')
\ No newline at end of file
diff --git a/tests/unit/main.py b/tests/unit/main.py
index df9783a..4cae687 100644
--- a/tests/unit/main.py
+++ b/tests/unit/main.py
@@ -7,6 +7,8 @@
 from test_dsmc.octree import *
 from test_dsmc.diagnostics import *
 from test_dsmc.common import *
+from test_dsmc.mesh.mesh2d import *
+from test_dsmc.boundary import *
 
 if __name__ == '__main__':
 	unittest.main()
diff --git a/tests/unit/test_dsmc/boundary.py b/tests/unit/test_dsmc/boundary.py
new file mode 100644
index 0000000..ed983fb
--- /dev/null
+++ b/tests/unit/test_dsmc/boundary.py
@@ -0,0 +1,195 @@
+import numpy as np
+import dsmc.boundary as bo
+import unittest
+
+class TestCommon(unittest.TestCase):
+    def test__check_if_parallel(self):
+        v1 = np.array([1.0, 0.0, 0.0])
+        v2 = np.array([1.0, 0.0, 0.0])
+        v3 = np.array([0.0, 1.0, 1.0])
+        v4 = np.array([1.0e-6, 10.0, 10.0])
+        
+        self.assertTrue(bo._check_if_parallel(v1, v2))
+        self.assertFalse(bo._check_if_parallel(v1, v3))
+        self.assertFalse(bo._check_if_parallel(v1, v4))
+        
+    def test__intersect(self):
+        l0 = np.array([-1.0, -1.0, 0.0])
+        l1 = np.array([1.0, 1.0, 0.0])
+        
+        p0 = np.array([0.0, -1.0, -1.0])
+        p1 = np.array([0.0, 1.0, -1.0])
+        p2 = np.array([0.0, -1.0, 1.0])
+        
+        intersected, n_l, n_p, t = bo._intersect(l0, l1, p0, p1, p2)
+        
+        self.assertTrue(intersected)
+        
+        self.assertEqual(2.0, n_l[0])
+        self.assertEqual(2.0, n_l[1])
+        self.assertEqual(0.0, n_l[2])
+        
+        self.assertEqual(4.0, n_p[0])
+        self.assertEqual(0.0, n_p[1])
+        self.assertEqual(0.0, n_p[2])
+        
+        self.assertEqual(0.5, t)
+        
+    def test__calc_nr(self):
+        l0 = np.array([-1.0, -1.0, 0.0])
+        l1 = np.array([1.0, 1.0, 0.0])
+        
+        p0 = np.array([0.0, -1.0, -1.0])
+        p1 = np.array([0.0, -1.0, 1.0])
+        p2 = np.array([0.0, 1.0, 1.0])
+        
+        intersected, n_l, n_p, t = bo._intersect(l0, l1, p0, p1, p2)
+        n_r = bo._calc_nr(n_l, n_p)
+        
+        self.assertEqual(-2.0, n_r[0])
+        self.assertEqual(2.0, n_r[1])
+        self.assertEqual(0.0, n_r[2])
+        
+    def test__reflect(self):
+        vel = np.array([123.0, 123.0, 0.0])
+        pos_old = np.array([-1.0, -1.0, 0.0])
+        pos = np.array([1.0, 1.0, 0.0])
+        
+        domain = np.array([(-2.0, 0.0), (-1.0, 1.0), (-1.0, 1.0)])
+        
+        p0, p1, p2 = bo._get_plane(domain, 0, 1)
+        
+        new_vel, new_pos, new_pos_old = bo._reflect(vel, pos, pos_old, p0, p1, p2, domain)
+        
+        self.assertEqual(-123.0, new_vel[0])
+        self.assertEqual(123.0, new_vel[1])
+        self.assertEqual(0.0, new_vel[2])
+        
+        self.assertEqual(-1.0, new_pos[0])
+        self.assertEqual(1.0, new_pos[1])
+        self.assertEqual(0.0, new_pos[2])
+        
+        self.assertEqual(0.0, new_pos_old[0])
+        self.assertEqual(0.0, new_pos_old[1])
+        self.assertEqual(0.0, new_pos_old[2])
+        
+    def test__get_plane1(self):
+        domain = np.array([(0, 1), (2, 4), (4, 7)])
+        axis = 0
+        
+        p0, p1, p2 = bo._get_plane(domain, axis, 0)
+        
+        self.assertEqual((3,), p0.shape)
+        self.assertEqual((3,), p1.shape)
+        self.assertEqual((3,), p2.shape)
+        
+        self.assertEqual(0.0, p0[0])
+        self.assertEqual(2.0, p0[1])
+        self.assertEqual(4.0, p0[2])
+        
+        self.assertEqual(0.0, p1[0])
+        self.assertEqual(2.0, p1[1])
+        self.assertEqual(7.0, p1[2])
+        
+        self.assertEqual(0.0, p2[0])
+        self.assertEqual(4.0, p2[1])
+        self.assertEqual(4.0, p2[2])
+        
+        p0, p1, p2 = bo._get_plane(domain, axis, 1)
+        
+        self.assertEqual((3,), p0.shape)
+        self.assertEqual((3,), p1.shape)
+        self.assertEqual((3,), p2.shape)
+        
+        self.assertEqual(1.0, p0[0])
+        self.assertEqual(2.0, p0[1])
+        self.assertEqual(4.0, p0[2])
+        
+        self.assertEqual(1.0, p1[0])
+        self.assertEqual(4.0, p1[1])
+        self.assertEqual(4.0, p1[2])
+        
+        self.assertEqual(1.0, p2[0])
+        self.assertEqual(2.0, p2[1])
+        self.assertEqual(7.0, p2[2])
+        
+    def test__get_plane2(self):
+        domain = np.array([(0, 1), (2, 4), (4, 7)])
+        axis = 1
+        
+        p0, p1, p2 = bo._get_plane(domain, axis, 0)
+        
+        self.assertEqual((3,), p0.shape)
+        self.assertEqual((3,), p1.shape)
+        self.assertEqual((3,), p2.shape)
+        
+        self.assertEqual(0.0, p0[0])
+        self.assertEqual(2.0, p0[1])
+        self.assertEqual(4.0, p0[2])
+        
+        self.assertEqual(1.0, p1[0])
+        self.assertEqual(2.0, p1[1])
+        self.assertEqual(4.0, p1[2])
+        
+        self.assertEqual(0.0, p2[0])
+        self.assertEqual(2.0, p2[1])
+        self.assertEqual(7.0, p2[2])
+        
+        p0, p1, p2 = bo._get_plane(domain, axis, 1)
+        
+        self.assertEqual((3,), p0.shape)
+        self.assertEqual((3,), p1.shape)
+        self.assertEqual((3,), p2.shape)
+        
+        self.assertEqual(0.0, p0[0])
+        self.assertEqual(4.0, p0[1])
+        self.assertEqual(4.0, p0[2])
+        
+        self.assertEqual(0.0, p1[0])
+        self.assertEqual(4.0, p1[1])
+        self.assertEqual(7.0, p1[2])
+        
+        self.assertEqual(1.0, p2[0])
+        self.assertEqual(4.0, p2[1])
+        self.assertEqual(4.0, p2[2])
+        
+    def test__get_plane3(self):
+        domain = np.array([(0, 1), (2, 4), (4, 7)])
+        axis = 2
+        
+        p0, p1, p2 = bo._get_plane(domain, axis, 0)
+        
+        self.assertEqual((3,), p0.shape)
+        self.assertEqual((3,), p1.shape)
+        self.assertEqual((3,), p2.shape)
+        
+        self.assertEqual(0.0, p0[0])
+        self.assertEqual(2.0, p0[1])
+        self.assertEqual(4.0, p0[2])
+        
+        self.assertEqual(0.0, p1[0])
+        self.assertEqual(4.0, p1[1])
+        self.assertEqual(4.0, p1[2])
+        
+        self.assertEqual(1.0, p2[0])
+        self.assertEqual(2.0, p2[1])
+        self.assertEqual(4.0, p2[2])
+        
+        p0, p1, p2 = bo._get_plane(domain, axis, 1)
+        
+        self.assertEqual((3,), p0.shape)
+        self.assertEqual((3,), p1.shape)
+        self.assertEqual((3,), p2.shape)
+        
+        self.assertEqual(0.0, p0[0])
+        self.assertEqual(2.0, p0[1])
+        self.assertEqual(7.0, p0[2])
+        
+        self.assertEqual(1.0, p1[0])
+        self.assertEqual(2.0, p1[1])
+        self.assertEqual(7.0, p1[2])
+        
+        self.assertEqual(0.0, p2[0])
+        self.assertEqual(4.0, p2[1])
+        self.assertEqual(7.0, p2[2])
+        
\ No newline at end of file
diff --git a/tests/unit/test_dsmc/mesh/__init__.py b/tests/unit/test_dsmc/mesh/__init__.py
new file mode 100644
index 0000000..633f866
--- /dev/null
+++ b/tests/unit/test_dsmc/mesh/__init__.py
@@ -0,0 +1,2 @@
+# -*- coding: utf-8 -*-
+
diff --git a/tests/unit/test_dsmc/mesh/mesh2d.py b/tests/unit/test_dsmc/mesh/mesh2d.py
new file mode 100644
index 0000000..e993b28
--- /dev/null
+++ b/tests/unit/test_dsmc/mesh/mesh2d.py
@@ -0,0 +1,143 @@
+import dsmc.mesh.mesh2d as msh
+import unittest
+import numpy as np
+
+class TestMesh2(unittest.TestCase):
+    def test_constructor(self):
+        mesh = msh.Mesh2d()
+        
+        self.assertEqual(1.0, mesh.cell_size1)
+        self.assertEqual(1.0, mesh.cell_size2)
+        self.assertEqual(0.0, mesh.min1)
+        self.assertEqual(0.0, mesh.min2)
+        self.assertEqual(1, mesh.n_cells1)
+        self.assertEqual(1, mesh.n_cells2)
+        self.assertEqual(msh.Plane.XY, mesh.plane)
+        
+    def test__get_cell_id1(self):
+        val1_f = -1.0
+        val1_t = 0.1
+        val2_f = -1.0
+        val2_t = 0.1
+        n_cells1 = 10
+        n_cells2 = 10
+        min1 = -0.5
+        min2 = -0.5
+        cell_size1 = 0.1
+        cell_size2 = 0.1
+        
+        res1_f = msh._get_cell_id(val1_t, val2_f, n_cells1, n_cells2, min1, min2, cell_size1, cell_size2)
+        res2_f = msh._get_cell_id(val1_f, val2_t, n_cells1, n_cells2, min1, min2, cell_size1, cell_size2)
+        res3_t = msh._get_cell_id(val1_t, val2_t, n_cells1, n_cells2, min1, min2, cell_size1, cell_size2)
+        
+        self.assertFalse(res1_f[0])
+        self.assertFalse(res2_f[0])
+        self.assertTrue(res3_t[0])
+        
+    def test__get_cell_id2(self):
+        val1 = (-0.45, -0.5)
+        val2 = (-0.35, -0.35)
+        n_cells1 = 10
+        n_cells2 = 10
+        min1 = -0.5
+        min2 = -0.5
+        cell_size1 = 0.1
+        cell_size2 = 0.1
+        
+        res1 = msh._get_cell_id(val1[0], val1[1], n_cells1, n_cells2, min1, min2, cell_size1, cell_size2)
+        res2 = msh._get_cell_id(val2[0], val2[1], n_cells1, n_cells2, min1, min2, cell_size1, cell_size2)
+        
+        self.assertTrue(res1[0])
+        self.assertTrue(res2[0])
+        
+        self.assertEqual(0, res1[1])
+        self.assertEqual(11, res2[1])
+        
+    def test_get_cell_id(self):
+        mesh = msh.Mesh2d()
+        
+        mesh.n_cells1 = 10
+        mesh.n_cells2 = 10
+        mesh.min1 = -0.5
+        mesh.min2 = -0.5
+        mesh.cell_size1 = 0.1
+        mesh.cell_size2 = 0.1
+        
+        pos = np.array([-0.45, -0.35, -0.25])
+        
+        # XY
+        mesh.plane = msh.Plane.XY
+        res = mesh.get_cell_id(pos)
+        
+        self.assertTrue(res[0])
+        self.assertEqual(10, res[1])
+        
+        # YZ
+        mesh.plane = msh.Plane.YZ
+        res = mesh.get_cell_id(pos)
+        
+        self.assertTrue(res[0])
+        self.assertEqual(21, res[1])
+        
+        # XZ
+        mesh.plane = msh.Plane.XZ
+        res = mesh.get_cell_id(pos)
+        
+        self.assertTrue(res[0])
+        self.assertEqual(20, res[1])
+        
+    def test__sort(self):
+        mesh = msh.Mesh2d()
+        
+        mesh.n_cells1 = 10
+        mesh.n_cells2 = 10
+        mesh.min1 = -0.5
+        mesh.min2 = -0.5
+        mesh.cell_size1 = 0.1
+        mesh.cell_size2 = 0.1
+        
+        pos1 = np.array([-0.45, -0.35, 0.0])
+        pos2 = np.array([-0.45, -3.50, 0.0])
+        positions = np.array([pos1, pos2])
+        
+        values1 = np.array([positions[i][0] for i in range(len(positions))])
+        values2 = np.array([positions[i][1] for i in range(len(positions))])
+        
+        inside, ids = msh._sort(values1, values2, mesh.n_cells1, mesh.n_cells2, mesh.min1, mesh.min2, mesh.cell_size1, mesh.cell_size2)
+        
+        self.assertTrue(inside[0])
+        self.assertFalse(inside[1])
+        
+        self.assertEqual(10, ids[0])
+        self.assertEqual(0, ids[1])
+        
+    def test_sort(self):
+        mesh = msh.Mesh2d()
+        
+        mesh.n_cells1 = 10
+        mesh.n_cells2 = 10
+        mesh.min1 = -0.5
+        mesh.min2 = -0.5
+        mesh.cell_size1 = 0.1
+        mesh.cell_size2 = 0.1
+        
+        positions = []
+        
+        positions.append([-0.45, -0.45, 0])
+        positions.append([-0.45, -0.45, 0])
+        
+        positions.append([0.45, 0.45, 0])
+        
+        mesh.sort(positions)
+        
+        self.assertEqual(len(mesh.cells), mesh.n_cells1 * mesh.n_cells2)
+        self.assertEqual(len(mesh.cells[0]), 2)
+        self.assertEqual(len(mesh.cells[99]), 1)
+        
+        self.assertEqual(mesh.cells[0][0], 0)
+        self.assertEqual(mesh.cells[0][1], 1)
+        self.assertEqual(mesh.cells[99][0], 2)
+        
+        for i in range(len(mesh.cells)):
+            if i != 0 and i != 99:
+                self.assertEqual(0, len(mesh.cells[i]))
\ No newline at end of file
diff --git a/tests/unit/test_dsmc/particles.py b/tests/unit/test_dsmc/particles.py
index 700c46f..0bed9aa 100644
--- a/tests/unit/test_dsmc/particles.py
+++ b/tests/unit/test_dsmc/particles.py
@@ -41,8 +41,9 @@ def test_calc_positions(self):
         x = (-1.0, 1.0)
         y = (2.0, 3.0)
         z = (-2.0, -1.0)
+        X = np.array([x, y, z])
         N = 1000
-        positions = pa.calc_positions(x, y, z, N)
+        positions = pa.calc_positions(X, N)
 
         self.assertEqual(N, len(positions))