diff --git a/map.py b/map.py
index 4adee05..fad2927 100644
--- a/map.py
+++ b/map.py
@@ -4,6 +4,8 @@ import numpy as np
 import sys
 import json
 
+from tqdm import tqdm, trange
+
 
 EPS = sys.float_info.epsilon
 
@@ -32,23 +34,30 @@ class Triangle:
         self.v3p = self.plane_coords(self.v3)
 
     def intersect(self, ray):
+        '''
+        :param ray: numpy.array(3) Ray cast from [0, 0, 0]
+        :return: A tuple of [bool, Optional[numpy.array(3)]]:
+                 The first component is True if the ray intersects the triangle, False otherwise.
+                 The second component marks the intersection of the triangle's plane (even if outside the triangle), or
+                 None if the ray does not intersect the plane at all.
+        '''
         # https://en.wikipedia.org/wiki/Möller–Trumbore_intersection_algorithm
         ray_cross_e2 = np.cross(ray, self.e2)
         dot = np.dot(self.e1, ray_cross_e2)
         if abs(dot) < EPS:
-            return None
+            return False, None
         s = -self.v1
         u = np.dot(s, ray_cross_e2) / dot
-        if (u < 0 and abs(u) > EPS) or (u > 1 and abs(u - 1) > EPS):
-            return None
         s_cross_e1 = np.cross(s, self.e1)
         v = np.dot(ray, s_cross_e1) / dot
-        if (v < 0 and abs(v) > EPS) or (u + v > 1 and abs(u + v - 1) > EPS):
-            return None
         t = np.dot(self.e2, s_cross_e1) / dot
+        if (u < 0 and abs(u) > EPS) or (u > 1 and abs(u - 1) > EPS):
+            return False, ray * t
+        if (v < 0 and abs(v) > EPS) or (u + v > 1 and abs(u + v - 1) > EPS):
+            return False, ray * t
         if t <= EPS:
-            return None
-        return ray * t
+            return False, None
+        return True, ray * t
 
     def plane_coords(self, v):
         u = np.dot(self.e1 / self.e1n, v - self.v1)
@@ -57,6 +66,16 @@ class Triangle:
         return np.array([u, v])
 
 
+class Face:
+
+    def __init__(self, vs):
+        self.vertices = vs
+        self.tris = []
+        for i in range(1, len(vs)-1):
+            self.tris.append(Triangle(vs[0], vs[i], vs[i+1]))
+        self.plane_coords = [self.tris[0].plane_coords(v) for v in self.vertices]
+
+
 def convert_to_cartesian(lon, lat):
     x = cos(lon*pi/180)
     y = sin(lon*pi/180)
@@ -65,52 +84,71 @@ def convert_to_cartesian(lon, lat):
     return np.array([x/mag, y/mag, z/mag])
 
 
-def map_poly(tri, poly):
+def map_poly(tri, poly, ref=None):
+    if ref is None:
+        ref = tri
     p0 = convert_to_cartesian(*poly[0])
-    p0i = tri.intersect(p0)
+    p0b, p0i = tri.intersect(p0)
     mapped = []
     new = True
     for i in range(1, len(poly)):
         p1 = convert_to_cartesian(*poly[i])
-        p1i = tri.intersect(p1)
-        if p0i is not None and p1i is not None:
+        p1b, p1i = tri.intersect(p1)
+        if p0b or p1b and p0i is not None and p1i is not None:
             if new:
-                mapped.append([tri.plane_coords(p0i), tri.plane_coords(p1i)])
+                mapped.append([ref.plane_coords(p0i), ref.plane_coords(p1i)])
                 new = False
             else:
-                mapped[-1].append(tri.plane_coords(p1i))
+                mapped[-1].append(ref.plane_coords(p1i))
         else:
             new = True
         p0 = p1
         p0i = p1i
+        p0b = p1b
     return mapped
 
 
 if __name__ == '__main__':
-    t = Triangle(np.array([5,0,0]), np.array([0,5,0]), np.array([0,0,5]))
-    countries = {}
+    fs = [Face([a, b, c, d]) for a, b, c, d in [
+          (np.array([1,1,1]), np.array([-1,1,1]), np.array([-1,-1,1]), np.array([1,-1,1])),  # top
+          (np.array([1,1,-1]), np.array([-1,1,-1]), np.array([-1,-1,-1]), np.array([1,-1,-1])),  # bottom
+          (np.array([-1,-1,-1]), np.array([-1,1,-1]), np.array([-1,1,1]), np.array([-1,-1,1])),  # left
+          (np.array([1,-1,-1]), np.array([1,1,-1]), np.array([1,1,1]), np.array([1,-1,1])),  # right
+          (np.array([-1,-1,-1]), np.array([1,-1,-1]), np.array([1,-1,1]), np.array([-1,-1,1])),  # front
+          (np.array([-1,1,-1]), np.array([1,1,-1]), np.array([1,1,1]), np.array([-1,1,1])),  # back
+    ]]
     with open('countries.geojson', 'r') as f:
         j = json.load(f)
-        for f in j['features']:
-            cc = f['properties']['ADMIN']
-            if f['geometry']['type'] == 'MultiPolygon':
-                countries[cc] = []
-                for poly in f['geometry']['coordinates']:
-                    countries[cc].extend(map_poly(t, poly[0]))
-            elif f['geometry']['type'] == 'Polygon':
-                countries[cc] = map_poly(tri, f['geometry']['coordinates'][0])
-    minx = min(t.v1p[0], t.v2p[0], t.v3p[0])
-    miny = min(t.v1p[1], t.v2p[1], t.v3p[1])
-    maxx = max(t.v1p[0], t.v2p[0], t.v3p[0])
-    maxy = max(t.v1p[1], t.v2p[1], t.v3p[1])
-    w = maxx - minx
-    h = maxy - miny
-    print(f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="{minx} {miny} {w} {h}" width="{w}" height="{h}">')
-    print(f'<path fill="none" stroke="black" stroke-width="0.01" d="M {t.v1p[0]} {t.v1p[1]} L {t.v2p[0]} {t.v2p[1]} L {t.v3p[0]} {t.v3p[1]} Z" />')
-    for c, l in countries.items():
-        for ls in l:
-            print(f'<path fill="none" stroke="red" stroke-width="0.01" d="M {ls[0][0]} {ls[0][1]}', end='')
-            for x, y in ls[1:]:
-                print(f' L {x} {y}', end='')
-            print('" />')
-    print('</svg>')
+
+    for i, face in tqdm(list(enumerate(fs)), desc='Faces   '):
+        countries = {}
+        for t in face.tris:
+            for f in tqdm(j['features'], desc='Features', total=len(face.tris)*len(j['features'])):
+                cc = f['properties']['ADMIN']
+                if f['geometry']['type'] == 'MultiPolygon':
+                    for poly in f['geometry']['coordinates']:
+                        countries.setdefault(cc, []).extend(map_poly(t, poly[0], face.tris[0]))
+                elif f['geometry']['type'] == 'Polygon':
+                    countries.setdefault(cc, []).extend(map_poly(t, f['geometry']['coordinates'][0], face.tris[0]))
+        minx = min(v[0] for v in face.plane_coords)
+        miny = min(v[1] for v in face.plane_coords)
+        maxx = max(v[0] for v in face.plane_coords)
+        maxy = max(v[1] for v in face.plane_coords)
+        w = maxx - minx
+        h = maxy - miny
+        with open(f'face{i}.svg', 'w') as svg:
+            svg.write(f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="{minx} {miny} {w} {h}" width="{w}" height="{h}" id="face{i}">\n')
+            svg.write(f'  <g id="cut">\n    <path id="outline" fill="none" stroke="black" stroke-width="0.01" d="M {face.plane_coords[0][0]} {face.plane_coords[0][1]}')
+            for v in face.plane_coords[1:]:
+                svg.write(f' L {v[0]} {v[1]}')
+            svg.write(' Z" />\n  </g>\n  <g id="engrave">\n')
+            for c, l in countries.items():
+                if l:
+                    svg.write(f'    <g id="{c}">\n')
+                    for ls in l:
+                        svg.write(f'      <path fill="none" stroke="red" stroke-width="0.01" d="M {ls[0][0]} {ls[0][1]}')
+                        for x, y in ls[1:]:
+                            svg.write(f' L {x} {y}')
+                        svg.write('" />\n')
+                    svg.write('    </g>\n')
+            svg.write('  </g>\n</svg>\n')