diff --git a/Cassiopee/XCore/XCore/common/mem.cpp b/Cassiopee/XCore/XCore/common/mem.cpp
index a52017271..84b522b28 100644
--- a/Cassiopee/XCore/XCore/common/mem.cpp
+++ b/Cassiopee/XCore/XCore/common/mem.cpp
@@ -75,3 +75,77 @@ void *xresize(void *ptr, E_Int nbytes, const char *file, E_Int line)
     }
     return ptr;
 }
+
+Mem_arena::Mem_arena()
+{
+    arena = NULL;
+    size = 0;
+    capacity = 0;
+    nresize = 0;
+    nalloc = 0;
+}
+
+Mem_arena::Mem_arena(size_t nbytes)
+{
+    arena = XMALLOC(nbytes);
+    size = 0;
+    capacity = nbytes;
+    nresize = 0;
+    nalloc = 0;
+}
+
+void Mem_arena::reserve(size_t nbytes)
+{
+    assert(arena == NULL);
+    arena = XMALLOC(nbytes);
+    capacity = nbytes;
+}
+
+void *Mem_arena::alloc(size_t nbytes, const char *file, int line)
+{
+    void *ptr = NULL;
+    
+    nalloc++;
+    
+    if (size + nbytes > capacity) {
+        /*
+        nresize++;
+        size_t new_size = size + nbytes;
+        size_t new_capacity = (size_t)((E_Float)new_size * 1.5);
+        arena = XRESIZE(arena, new_capacity);
+        assert(size + nbytes < new_capacity);
+        
+        ptr = (void *)((char *)arena + size);
+        
+        size = new_size;
+        capacity = new_capacity;
+        */
+        fprintf(stderr, "Mem_arena: out of memory at %s:%d\n", file, line);
+        print_stats();
+        exit(1);
+    } else {
+        ptr = (void *)((char *)arena + size);
+        size += nbytes;
+    }
+    
+    return ptr;
+}
+
+void Mem_arena::drop()
+{
+    XFREE(arena);
+    size = capacity = nalloc = nresize = 0;
+}
+
+void Mem_arena::print_stats()
+{
+    printf("\n");
+    printf("    /**** Mem_arena stats ****/\n");
+    printf("      size:     %lu B\n", size);
+    printf("      capacity: %lu B\n", capacity);
+    printf("      usage:    %.2f %%\n", size * 100.0f / capacity);
+    printf("      number of allocations:   %zu\n", nalloc);
+    printf("      number of reallocations: %zu\n", nresize);
+    printf("    /*************************/\n");
+    printf("\n");
+}
\ No newline at end of file
diff --git a/Cassiopee/XCore/XCore/common/mem.h b/Cassiopee/XCore/XCore/common/mem.h
index 992af2fa4..16c202fe8 100644
--- a/Cassiopee/XCore/XCore/common/mem.h
+++ b/Cassiopee/XCore/XCore/common/mem.h
@@ -39,4 +39,28 @@ void *xcalloc(E_Int, E_Int, const char *, E_Int);
 void xfree(void *, const char *, E_Int);
 void *xresize(void *, E_Int, const char *, E_Int);
 
+struct Mem_arena {
+    void *arena;
+
+    size_t size;
+    size_t capacity;
+
+    size_t nresize;
+    size_t nalloc;
+
+    Mem_arena();
+
+    Mem_arena(size_t nbytes);
+
+    void reserve(size_t nbytes);
+
+    void *alloc(size_t nbytes, const char *file, int line);
+
+    void drop();
+
+    void print_stats();
+};
+
+#define ARENA_ALLOC(nbytes) arena.alloc((nbytes), __FILE__, __LINE__)
+
 #endif
diff --git a/Cassiopee/XCore/XCore/intersectMesh/AABB.cpp b/Cassiopee/XCore/XCore/intersectMesh/AABB.cpp
new file mode 100644
index 000000000..89bc83981
--- /dev/null
+++ b/Cassiopee/XCore/XCore/intersectMesh/AABB.cpp
@@ -0,0 +1,31 @@
+#include "AABB.h"
+#include "mesh.h"
+
+AABB::AABB()
+: xmin(EFLOATMAX), ymin(EFLOATMAX), zmin(EFLOATMAX),
+  xmax(EFLOATMIN), ymax(EFLOATMIN), zmax(EFLOATMIN)
+{}
+
+AABB::AABB(const IMesh &M, E_Int *ids, E_Int count)
+{
+    xmin = ymin = zmin = EFLOATMAX;
+    xmax = ymax = zmax = EFLOATMIN;
+
+    for (E_Int i = 0; i < count; i++) {
+        E_Int fid = ids[i];
+
+        const auto &pn = M.F[fid];
+
+        for (E_Int p : pn) {
+            E_Float x = M.X[p];
+            E_Float y = M.Y[p];
+            E_Float z = M.Z[p];
+            if (x < xmin) xmin = x;
+            if (y < ymin) ymin = y;
+            if (z < zmin) zmin = z;
+            if (x > xmax) xmax = x;
+            if (y > ymax) ymax = y;
+            if (z > zmax) zmax = z;
+        }
+    }
+}
\ No newline at end of file
diff --git a/Cassiopee/XCore/XCore/intersectMesh/AABB.h b/Cassiopee/XCore/XCore/intersectMesh/AABB.h
index e760329c2..a07c133d0 100644
--- a/Cassiopee/XCore/XCore/intersectMesh/AABB.h
+++ b/Cassiopee/XCore/XCore/intersectMesh/AABB.h
@@ -1,10 +1,18 @@
 #pragma once
 
+#include "xcore.h"
+
+struct IMesh;
+
 struct AABB {
-    E_Float xmin = EFLOATMAX;
-    E_Float xmax = EFLOATMIN;
-    E_Float ymin = EFLOATMAX;
-    E_Float ymax = EFLOATMIN;
-    E_Float zmin = EFLOATMAX;
-    E_Float zmax = EFLOATMIN;
+    E_Float xmin;
+    E_Float ymin;
+    E_Float zmin;
+    E_Float xmax;
+    E_Float ymax;
+    E_Float zmax;
+
+    AABB();
+
+    AABB(const IMesh &M, E_Int *ids, E_Int count);
 };
\ No newline at end of file
diff --git a/Cassiopee/XCore/XCore/intersectMesh/BVH.cpp b/Cassiopee/XCore/XCore/intersectMesh/BVH.cpp
new file mode 100644
index 000000000..56dd3b545
--- /dev/null
+++ b/Cassiopee/XCore/XCore/intersectMesh/BVH.cpp
@@ -0,0 +1,107 @@
+#include "BVH.h"
+#include "mesh.h"
+
+// Creates the tree
+
+static size_t leaf_count = 0;
+
+BVH_node *BVH_create_node(AABB aabb, BVH_node *left_child, BVH_node *right_child,
+    E_Int *ids, E_Int count, Mem_arena &arena)
+{
+    BVH_node *node = (BVH_node *)ARENA_ALLOC(sizeof(BVH_node));
+
+    node->bbox = aabb;
+    node->left = left_child;
+    node->right = right_child;
+    node->elements = (E_Int *)ARENA_ALLOC(count * sizeof(E_Int));
+    assert(count > 0);
+    memcpy(node->elements, ids, count * sizeof(E_Int));
+
+    leaf_count++;
+
+    return node;
+}
+
+
+BVH_node *BVH_create_node(AABB aabb, BVH_node *left_child,
+    BVH_node *right_child, Mem_arena &arena)
+{
+    BVH_node *node = (BVH_node *)ARENA_ALLOC(sizeof(BVH_node));
+
+    node->bbox = aabb;
+    node->left = left_child;
+    node->right = right_child;
+    node->elements = NULL;
+
+    return node;
+}
+
+BVH_node *BVH_create_node(E_Int *ids, E_Int count, const IMesh &M,
+    Mem_arena &arena)
+{
+    AABB aabb(M, ids, count);
+
+    if (count <= MAX_FACES_PER_LEAF) {
+        return BVH_create_node(aabb, NULL, NULL, ids, count, arena);
+    }
+
+    E_Float dx = aabb.xmax - aabb.xmin;
+    E_Float dy = aabb.ymax - aabb.ymin;
+    E_Float dz = aabb.zmax - aabb.zmin;
+
+    E_Float *pCoor;
+
+    if (dx >= dy && dx >= dz) {
+        pCoor = (E_Float *)M.X.data();
+    } else if (dy >= dz) {
+        assert(dy >= dx);
+        pCoor = (E_Float *)M.Y.data();
+    } else {
+        assert(dz >= dx && dz >= dy);
+        pCoor = (E_Float *)M.Z.data();
+    }
+
+    std::sort(ids, ids + count, [&] (E_Int i, E_Int j) 
+    {
+        E_Float cci = 0.0, ccj = 0.0;
+        
+        const auto &pni = M.F[i];
+        for (const E_Int p : pni) cci += pCoor[p];
+        cci /= pni.size();
+
+        const auto &pnj = M.F[j];
+        for (const E_Int p : pnj) ccj += pCoor[p];
+        ccj /= pnj.size();
+        
+        return cci < ccj;
+    });
+
+    BVH_node *left = BVH_create_node(ids, count/2, M, arena);
+    BVH_node *right = BVH_create_node(ids + count/2, count - count/2, M, arena);
+
+    return BVH_create_node(aabb, left, right, arena);
+}
+
+BVH::BVH(const IMesh &M)
+{
+    if (M.skin.empty()) {
+        puts("empty skin!");
+        exit(1);
+    }
+
+    size_t NF = M.skin.size();
+
+    printf("Number of faces: %zu\n", NF);
+
+    arena.reserve(10 * NF * sizeof(E_Int));
+
+    E_Int *ids = (E_Int *)ARENA_ALLOC(NF * sizeof(E_Int));
+
+    for (size_t i = 0; i < NF; i++) ids[i] = M.skin[i];
+
+    root = BVH_create_node(ids, NF, M, arena);
+
+    printf("Leaves: %zu\n", leaf_count);
+
+    arena.print_stats();
+}
\ No newline at end of file
diff --git a/Cassiopee/XCore/XCore/intersectMesh/BVH.h b/Cassiopee/XCore/XCore/intersectMesh/BVH.h
new file mode 100644
index 000000000..ddcd60023
--- /dev/null
+++ b/Cassiopee/XCore/XCore/intersectMesh/BVH.h
@@ -0,0 +1,32 @@
+#pragma once
+
+#include "AABB.h"
+#include "common/mem.h"
+
+#define MAX_FACES_PER_LEAF 10
+
+struct IMesh;
+
+struct BVH_node {
+    AABB bbox;
+    BVH_node *left;
+    BVH_node *right;
+    E_Int *elements;
+};
+
+BVH_node *BVH_create_node(AABB aabb, BVH_node *left_child,
+    BVH_node *right_child, Mem_arena &arena);
+
+BVH_node *BVH_create_node(AABB aabb, BVH_node *left_child, BVH_node *right_child,
+        E_Int *ids, E_Int count, Mem_arena &arena);
+
+BVH_node *BVH_create_node(E_Int *ids, E_Int count, const IMesh &M,
+    Mem_arena &arena);
+
+struct BVH {
+    BVH_node *root;
+
+    Mem_arena arena;
+
+    BVH(const IMesh &M);
+};
\ No newline at end of file
diff --git a/Cassiopee/XCore/XCore/intersectMesh/IntersectMesh_Init.cpp b/Cassiopee/XCore/XCore/intersectMesh/IntersectMesh_Init.cpp
index e242ade05..ca555d7cf 100644
--- a/Cassiopee/XCore/XCore/intersectMesh/IntersectMesh_Init.cpp
+++ b/Cassiopee/XCore/XCore/intersectMesh/IntersectMesh_Init.cpp
@@ -18,6 +18,7 @@
 */
 #include "mesh.h"
 #include "karray.h"
+//#include "BVH.h"
 
 PyObject *K_XCORE::IntersectMesh_Init(PyObject *self, PyObject *args)
 {
@@ -39,6 +40,13 @@ PyObject *K_XCORE::IntersectMesh_Init(PyObject *self, PyObject *args)
     // Init mesh
 
     IMesh *M = new IMesh(*karray.cn, karray.X, karray.Y, karray.Z, karray.npts);
+    /*
+    M->make_skin();
+    puts("Making BVH");
+    BVH bvh(*M);
+    puts("Done");
+    bvh.arena.drop();
+    */
 
     if (TAGS != Py_None) {
         E_Int size, nfld;
@@ -67,4 +75,4 @@ PyObject *K_XCORE::IntersectMesh_Init(PyObject *self, PyObject *args)
     PyObject *hook = PyCapsule_New((void *)M, "IntersectMesh", NULL);
 
     return hook;
-}
\ No newline at end of file
+}
diff --git a/Cassiopee/XCore/XCore/intersectMesh/dcel.cpp b/Cassiopee/XCore/XCore/intersectMesh/dcel.cpp
index d777f8980..01d94bcd0 100644
--- a/Cassiopee/XCore/XCore/intersectMesh/dcel.cpp
+++ b/Cassiopee/XCore/XCore/intersectMesh/dcel.cpp
@@ -1040,7 +1040,7 @@ void Dcel::trace_hedge(Hedge *sh, const Smesh &M, const Smesh &S, E_Int hid)
 
                     } else {
                         
-                        // E_Intersection, move
+                        // Intersection, move
 
                         current_face = h->twin->left;
     
diff --git a/Cassiopee/XCore/srcs.py b/Cassiopee/XCore/srcs.py
index e0124bcff..63a8a5686 100644
--- a/Cassiopee/XCore/srcs.py
+++ b/Cassiopee/XCore/srcs.py
@@ -20,6 +20,9 @@
             'XCore/common/mem.cpp',
             'XCore/common/common.cpp',
 
+            'XCore/intersectMesh/BVH.cpp',
+            'XCore/intersectMesh/AABB.cpp',
+
             'XCore/intersectMesh/extract.cpp',
 
             'XCore/intersectMesh/IntersectMesh_Init.cpp',