diff --git a/web/app/3d/debug.js b/web/app/3d/debug.js
index 2c1a9889..39dbb83d 100644
--- a/web/app/3d/debug.js
+++ b/web/app/3d/debug.js
@@ -64,11 +64,10 @@ function addGlobalDebugActions(app) {
       app.viewer.render();
     },
     AddHalfEdge: (he, color) => {
-      const points = [he.vertexA.point];
-      if (he.edge && he.edge.curve) {
-        he.edge.curve.approximate(10, he.vertexA.point, he.vertexB.point, points);
+      const points = he.edge.curve.verb.tessellate().map(p => new Vector().set3(p));
+      if (he.inverted) {
+        points.reverse();
       }
-      points.push(he.vertexB.point);
       window.__DEBUG__.AddPolyLine(points, color);  
     },
     AddFace: (face, color) => {
diff --git a/web/app/3d/modeler-app.js b/web/app/3d/modeler-app.js
index d0cc52f5..64ecad2e 100644
--- a/web/app/3d/modeler-app.js
+++ b/web/app/3d/modeler-app.js
@@ -88,11 +88,11 @@ App.prototype.scratchCode = function() {
   const box1 = app.TPI.brep.primitives.box(500, 500, 500);
   const box2 = app.TPI.brep.primitives.box(250, 250, 750, new Matrix3().translate(25, 25, 0));
   const box3 = app.TPI.brep.primitives.box(150, 600, 350, new Matrix3().translate(25, 25, -250));
-  let result = app.TPI.brep.bool.union(box1, box2);
-  // let result = app.TPI.brep.bool.subtract(box1, box2);
-  // result = app.TPI.brep.bool.subtract(result, box3);
-  app.addShellOnScene(box1);
-  app.addShellOnScene(box2);
+  // let result = app.TPI.brep.bool.union(box1, box2);
+  let result = app.TPI.brep.bool.subtract(box1, box2);
+  result = app.TPI.brep.bool.subtract(result, box3);
+  // app.addShellOnScene(box1);
+  app.addShellOnScene(result);
   app.viewer.render();
 };
 
diff --git a/web/app/3d/scene/brep-scene-object.js b/web/app/3d/scene/brep-scene-object.js
index 08beaf8a..359238ac 100644
--- a/web/app/3d/scene/brep-scene-object.js
+++ b/web/app/3d/scene/brep-scene-object.js
@@ -41,24 +41,18 @@ export class BREPSceneSolid extends SceneSolid {
 
   createEdges() {
     const visited = new Set();
-    for (let face of this.shell.faces) {
-      for (let halfEdge of face.outerLoop.halfEdges) {
-        if (!visited.has(halfEdge.edge)) {
-          visited.add(halfEdge.edge);
-          if (halfEdge.edge.data[EDGE_AUX] === undefined) {
-            const line = new THREE.Line(undefined, WIREFRAME_MATERIAL);
-            const contour = [halfEdge.vertexA.point];
-            halfEdge.edge.curve.approximate(10, halfEdge.vertexA.point, halfEdge.vertexB.point, contour);
-            contour.push(halfEdge.vertexB.point);
-            for (let p of contour) {
-              line.geometry.vertices.push(p.three());
-            }
-            this.wireframeGroup.add(line);
-            line.__TCAD_EDGE = halfEdge.edge;
-            halfEdge.edge.data['scene.edge'] = line;
-          }
+    for (let edge of this.shell.edges) {
+      if (edge.data[EDGE_AUX] === undefined) {
+        const line = new THREE.Line(undefined, WIREFRAME_MATERIAL);
+        const contour = edge.curve.verb.tessellate();
+        for (let p of contour) {
+          line.geometry.vertices.push(new THREE.Vector3().fromArray(p));
         }
+        this.wireframeGroup.add(line);
+        line.__TCAD_EDGE = edge;
+        edge.data['scene.edge'] = line;
       }
+
     }
   }
 
diff --git a/web/app/3d/tess/brep-tess.js b/web/app/3d/tess/brep-tess.js
index 1bdfd4a5..7a9fee4f 100644
--- a/web/app/3d/tess/brep-tess.js
+++ b/web/app/3d/tess/brep-tess.js
@@ -55,13 +55,7 @@ export default function(face) {
   }
   analyzeCurvature(face.surface.verb, triangles);
 
-  if (face.surface.inverted) {
-    triangles.forEach(t => t.reverse());
-  }
-
   return triangles.map(t => t.map(p => face.surface.point(p[0], p[1])));
-
-
 }
 
 function analyzeCurvature(nurbs, triangles) {
diff --git a/web/app/brep/brep-builder.js b/web/app/brep/brep-builder.js
index ae9caf43..725520ce 100644
--- a/web/app/brep/brep-builder.js
+++ b/web/app/brep/brep-builder.js
@@ -112,193 +112,6 @@ function assemble(walls, basePlane, lidPlane) {
   return shell;
 }
 
-function assembleRevolved(walls, baseSurface, lidSurface) {
-  const baseLoop = new Loop();
-  const lidLoop = new Loop();
-  const shell = new Shell();
-
-  // walls.reverse();
-  iterateSegments(walls, (wall, next) => {
-
-    const nullFace = new Face();
-    nullFace.outerLoop = new Loop();
-    function addEdge(he) {
-      he.edge = new Edge(Line.fromSegment(he.vertexA.point, he.vertexB.point))
-      const twin = new HalfEdge().setAB(he.vertexB, he.vertexA);
-      twin.loop = nullFace.outerLoop;
-      nullFace.outerLoop.halfEdges.push(twin);
-      he.edge.link(twin, he);
-    }
-
-    for (let he of wall.topEdges) {
-      addEdge(he);
-//      __DEBUG__.AddHalfEdge(he);
-    }
-
-    for (let i = next.bottomEdges.length - 1; i >= 0; i--) {
-      let he = next.bottomEdges[i];
-      addEdge(he);
-//      __DEBUG__.AddHalfEdge(he, 0xffffff);
-    }
-
-
-    linkSegments(nullFace.outerLoop.halfEdges);
-    nullFace.data.NULL_FACE = {
-      curve: vertIsoCurve(wall.surface, 1, true),
-      start: wall.topEdges[0].twin()
-    };
-//    __DEBUG__.AddPoint.apply(null, nullFace.data.NULL_FACE.curve.point(0.0))
-//    __DEBUG__.AddPoint.apply(null, nullFace.data.NULL_FACE.curve.point(0.1))
-//    __DEBUG__.AddPoint.apply(null, nullFace.data.NULL_FACE.curve.point(0.2))
-//    __DEBUG__.AddPoint.apply(null, nullFace.data.NULL_FACE.curve.point(0.3))
-//    __DEBUG__.AddPoint.apply(null, nullFace.data.NULL_FACE.curve.point(0.4))
-//    __DEBUG__.AddPoint.apply(null, nullFace.data.NULL_FACE.curve.point(0.5))
-    mergeNullFace(nullFace.data.NULL_FACE);
-  });
-
-
-  for (let wall of walls) {
-
-    function addHalfEdges(loop, edges) {
-//      for (let i = edges.length - 1; i >= 0; i--) {
-//        let he = edges[i];
-      for (let he of edges) {
-        he.edge = new Edge(Line.fromSegment(he.vertexA.point, he.vertexB.point))
-        const twin = new HalfEdge().setAB(he.vertexB, he.vertexA);
-        __DEBUG__.AddHalfEdge(twin)
-        twin.loop = loop;
-        loop.halfEdges.push(twin);
-        he.edge.link(twin, he);
-      }
-    }
-    addHalfEdges(baseLoop, wall.leftEdges);
-    addHalfEdges(lidLoop, wall.rightEdges);
-
-    for (let wallFace of wall.faces) {
-      shell.faces.push(wallFace);
-    }
-  }
-  
-  // walls.reverse();
-
-  lidLoop.halfEdges.reverse();
-  linkSegments(baseLoop.halfEdges);
-  linkSegments(lidLoop.halfEdges);
-
-  const baseFace = createFace(baseSurface, baseLoop);
-  const lidFace = createFace(lidSurface, lidLoop);
-
-  shell.faces.push(baseFace, lidFace);
-  shell.faces.forEach(f => f.shell = shell);
-  return shell;
-}
-
-function connectWalls(walls) {
-
-  iterateSegments(walls, (a, b) => {
-
-    function connect(halfEdgeA, halfEdgeB) {
-      const curve = Line.fromSegment(halfEdgeA.vertexA.point, halfEdgeA.vertexB.point);
-      new Edge(curve).link(halfEdgeA, halfEdgeB);
-    }
-
-    let aEdges = a.leftEdges;
-    let bEdges = b.rightEdges;
-
-    if (aEdges.length == 1 && bEdges.length == 1) {
-      connect(aEdges[0], bEdges[0])
-    } else {
-      throw "unsupported: use 'null-face' like it's done for the revolve and then merge it";
-    }
-  });
-}
-
-export function revolveToWallNurbs(basePath, surface, p0, v, angle) {
-  const nurbses = [];
-  const n = basePath.points.length;
-  for (let i = 0; i < n; i++) {
-    const curve = basePath.groups[i].toNurbs(surface);
-    const nurbs = new verb.geom.RevolvedSurface(curve.verb, p0.data(), v.data(), -angle);
-    nurbses.push(nurbs);
-  }
-  return nurbses;
-}
-
-function swap(obj, prop1, prop2) {
-  const tmp = obj[prop1];
-  obj[prop1] = obj[prop2];
-  obj[prop2] = tmp;
-}
-
-export function revolve(basePath, baseSurface, p0, v, angle) {
-
-  angle = -angle;
-
-  const baseLoop = new Loop();
-
-  const shell = new Shell();
-  const walls = [];
-
-  const n = basePath.points.length;
-
-  const baseVertices = [];
-  const lidVertices = [];
-
-  const nurbses = revolveToWallNurbs(basePath, baseSurface, p0, v, -angle);
-
-  for (let nurbs of nurbses) {
-    const domU = nurbs.domainU();
-    const domV = nurbs.domainV();
-    // profile of revolving becomes V direction
-    baseVertices.push(new Vertex(new Point().set3(nurbs.point(domU.min, domV.min))));
-    lidVertices.push(new Vertex(new Point().set3(nurbs.point(domU.max, domV.min))));
-  }
-
-  for (let i = 0; i < n; i++) {
-    let j = (i + 1) % n;
-    const nurbs = nurbses[i];
-    const wall = wallFromNUBRS(nurbs, false, baseVertices[i], lidVertices[i], lidVertices[j], baseVertices[j]);
-    walls.push(wall);
-  }
-
-  const normal = cad_utils.normalOfCCWSeq([lidVertices[2].point, lidVertices[1].point, lidVertices[0].point]);
-  const w = lidVertices[0].point.dot(normal);
-  const planeLid = new Plane(normal, w);
-
-  let revolved = assembleRevolved(walls, baseSurface, planeLid);
-  if (angle < 0) {
-    invert(revolved);
-  }
-  return revolved;
-}
-
-function createTwin(halfEdge) {
-  const twin = new HalfEdge();
-  twin.vertexA = halfEdge.vertexB;
-  twin.vertexB = halfEdge.vertexA;
-  twin.edge = halfEdge.edge;
-  if (halfEdge.edge.halfEdge1 == halfEdge) {
-    halfEdge.edge.halfEdge2 = twin;
-  }  else {
-    halfEdge.edge.halfEdge1 = twin;
-  }
-  return twin;
-}
-
-function createFace(surface, loop) {
-  const face = new Face(surface);
-  face.outerLoop = loop;
-  loop.face = face;
-  return face;
-}
-
-
-function createPlaneForLoop(normal, loop) {
-  const w = loop.halfEdges[0].vertexA.point.dot(normal);
-  const plane = new Plane(normal, w);
-  return plane;
-}
-
 function createBoundingNurbs(points, plane) {
   if (!plane) {
     const normal = cad_utils.normalOfCCWSeq(points);
@@ -322,55 +135,6 @@ function createBoundingNurbs(points, plane) {
   return nurbs;
 }
 
-export function linkHalfEdges(edge, halfEdge1, halfEdge2) {
-  halfEdge1.edge = edge;
-  halfEdge2.edge = edge;
-  edge.halfEdge1 = halfEdge1;
-  edge.halfEdge2 = halfEdge2;
-}
-
-export function createHalfEdge(loop, vertexA, vertexB) {
-  const halfEdge = new HalfEdge();
-  halfEdge.loop = loop;
-  halfEdge.vertexA = vertexA;
-  halfEdge.vertexB = vertexB;
-  loop.halfEdges.push(halfEdge);
-  return halfEdge;
-}
-
-export function linkSegments(halfEdges) {
-  iterateSegments(halfEdges, (prev, next) => {
-    prev.next = next;
-    next.prev = prev;
-  });
-}
-
-export function point(x, y, z) {
-  return new Point(x, y, z);
-}
-
-export function iterateSegments(items, callback) {
-  let length = items.length;
-  for (let i = 0; i < length; i++) {
-    let j = (i + 1) % length;
-    callback(items[i], items[j], i, j);
-  }
-}
-
-export function createPlaneLoop(vertices, curves) {
-
-  const loop = new Loop();
-
-  iterateSegments(vertices, (a, b, i) => {
-    const halfEdge = createHalfEdge(loop, a, b);
-    halfEdge.edge = new Edge(curves[i] ? curves[i] : Line.fromSegment(a.point, b.point));
-    return halfEdge;
-  });
-
-  linkSegments(loop.halfEdges);
-  return loop;
-}
-
 function bothClassOf(o1, o2, className) {
   return o1.constructor.name === className && o2.constructor.name === className;
 }
diff --git a/web/app/brep/brep-primitives.js b/web/app/brep/brep-primitives.js
index 46dda57f..092c11ac 100644
--- a/web/app/brep/brep-primitives.js
+++ b/web/app/brep/brep-primitives.js
@@ -1,4 +1,5 @@
-import * as BREPBuilder from './brep-builder'
+import {Point} from './geom/point'
+import {createPrism} from './brep-builder'
 import {Matrix3} from '../math/l3space'
 
 export function box(w, h, d, tr) {
@@ -8,11 +9,11 @@ export function box(w, h, d, tr) {
   if (!tr) {
     tr = IDENTITY;
   }
-  return BREPBuilder.createPrism([
-    tr._apply(BREPBuilder.point(-wh, -hh, dh)),
-    tr._apply(BREPBuilder.point( wh, -hh, dh)),
-    tr._apply(BREPBuilder.point( wh,  hh, dh)),
-    tr._apply(BREPBuilder.point(-wh,  hh, dh))
+  return createPrism([
+    tr._apply(new Point(-wh, -hh, dh)),
+    tr._apply(new Point( wh, -hh, dh)),
+    tr._apply(new Point( wh,  hh, dh)),
+    tr._apply(new Point(-wh,  hh, dh))
   ], d);
 }
 
diff --git a/web/app/brep/geom/impl/nurbs.js b/web/app/brep/geom/impl/nurbs.js
index 66506736..2708e9ab 100644
--- a/web/app/brep/geom/impl/nurbs.js
+++ b/web/app/brep/geom/impl/nurbs.js
@@ -17,45 +17,12 @@ export class NurbsCurve extends Curve {
     return new NurbsCurve(this.verb.transform(tr));
   }
   
-  approximate(resolution, from, to, out) {
-    const chunks = this.verb.divideByArcLength(10);
-    let startU = this.verb.closestParam(from.toArray());
-    let endU = this.verb.closestParam(to.toArray());
-    const reverse = startU > endU;
-    if (reverse) {
-      const tmp = startU;
-      startU = endU;
-      endU = tmp;
-      chunks.reverse();
-    }
-
-    for (let sample of chunks) {
-      const u = sample.u;
-      if (u > startU + math.TOLERANCE && u < endU - math.TOLERANCE) {
-        out.push(new Point().set3(this.verb.point(u)));
-      }
-    }
-  }
-
-  approximateU(resolution, paramFrom, paramTo, consumer) {
-    let u = paramFrom;
-    let endU = paramTo;
-    let step = this.verb.paramAtLength(resolution);
-    if (u > endU) {
-      step *= -1;
-    }
-    u += step;
-    for (;step > 0 ? u < endU : u > endU; u += step) {
-      consumer(u);
-    }
-  }
-  
   tangentAtPoint(point) {
-    return new Point().set3(this.verb.tangent(this.verb.closestParam(point.data())))._normalize();
+    return pt(this.verb.tangent(this.verb.closestParam(point.data())))._normalize();
   }
 
   tangentAtParam(param) {
-    return new Point().set3(this.verb.tangent(param ))._normalize();
+    return pt(this.verb.tangent(param ))._normalize();
   }
   
   closestDistanceToPoint(point) {
@@ -64,7 +31,7 @@ export class NurbsCurve extends Curve {
   }
 
   split(point) {
-    return this.verb.split(this.verb.closestParam(point.data)).map(v => new NurbsCurve(v));
+    return this.verb.split(this.verb.closestParam(point.data())).map(v => new NurbsCurve(v));
   }
 
   invert() {
@@ -72,17 +39,55 @@ export class NurbsCurve extends Curve {
   }
   
   point(u) {
-    return new Point().set3(this.verb.point(u));
+    return pt(this.verb.point(u));
   }
 
   intersectCurve(other, tol) {
-    return verb.geom.Intersect.curves(this.verb, other.verb, tol).map( i => ({
+    let isecs = [];
+    tol = tol || 1e-6;
+
+    const eq = (v1, v2) => math.areVectorsEqual3(v1, v2, tol);
+
+    function add(i0) {
+      for (let i1 of isecs) {
+        if (eq(i0.p0, i1.p0)) {
+          return;    
+        }    
+      }  
+      isecs.push(i0);
+    }
+
+    function isecOn(c0, c1, u0) {
+      const p0 = c0.verb.point(u0);
+      const u1 = c1.verb.closestParam(p0);
+      const p1 = c1.verb.point(u1);
+      if (eq(p0, p1)) {
+        if (c0 === other) {
+          add({u0: u1, u1: u0, p0: p1, p1: p0});
+        } else {
+          add({u0, u1, p0, p1});
+        }
+       
+      }
+    }
+    
+    isecOn(this, other, 0);
+    isecOn(this, other, 1);
+    isecOn(other, this, 0);
+    isecOn(other, this, 1);
+
+    verb.geom.Intersect.curves(this.verb, other.verb, tol).forEach( i => add({
       u0: i.u0,
       u1: i.u1,
-      p0: new Vector().set3(i.point0),
-      p1: new Vector().set3(i.point1)
+      p0: i.point0,
+      p1: i.point1
     }));
-  }
+    isecs.forEach(i => {
+      i.p0 = pt(i.p0);
+      i.p1 = pt(i.p1);
+    })
+    return isecs;
+}
 
   static createByPoints(points, degeree) {
     points = points.map(p => p.data());
@@ -112,7 +117,7 @@ export class NurbsSurface extends Surface {
 
   normal(point) {
     let uv = this.verb.closestParam(point.data());
-    let normal = new Vector().set3(this.verb.normal(uv[0], uv[1]));
+    let normal = pt(this.verb.normal(uv[0], uv[1]));
     if (this.inverted) {
       normal._negate();
     }
@@ -121,7 +126,7 @@ export class NurbsSurface extends Surface {
   }
 
   normalUV(u, v) {
-    let normal = new Vector().set3(this.verb.normal(u, v));
+    let normal = pt(this.verb.normal(u, v));
     if (this.inverted) {
       normal._negate();
     }
@@ -134,11 +139,11 @@ export class NurbsSurface extends Surface {
   }
 
   point(u, v) {
-    return new Point().set3(this.verb.point(u, v));
+    return pt(this.verb.point(u, v));
   }
 
   intersectSurfaceForSameClass(other, tol) {
-    const curves = verb.geom.Intersect.surfaces(this.verb, other.verb, tol);
+    const curves = verb_isec(this.verb, other.verb);
     let inverted = this.inverted !== other.inverted;
     return curves.map(curve => new NurbsCurve(inverted ?  curve.reverse() : curve));
   }
@@ -162,4 +167,32 @@ export class NurbsSurface extends Surface {
   isoCurveAlignV(param) {
     return this.isoCurve(param, false);
   }
-}
\ No newline at end of file
+}
+
+function dist(p1, p2) {
+  return math.distance3(p1[0], p1[1], p1[2], p2[0], p2[1], p2[2]);
+}
+
+function pt(data) {
+  return new Point().set3(data);
+}
+
+
+function verb_isec(nurbs1, nurbs2) {
+  const tol = 1e-3
+  const surface0 = nurbs1.asNurbs();
+  const surface1 = nurbs2.asNurbs();
+	var tess1 = verb.eval.Tess.rationalSurfaceAdaptive(surface0);
+	var tess2 = verb.eval.Tess.rationalSurfaceAdaptive(surface1);
+	var resApprox = verb.eval.Intersect.meshes(tess1,tess2);
+	var exactPls = resApprox.map(function(pl) {
+		return pl.map(function(inter) {
+			return verb.eval.Intersect.surfacesAtPointWithEstimate(surface0,surface1,inter.uv0,inter.uv1,tol);
+		});
+	});
+	return exactPls.map(function(x) {
+		return verb.eval.Make.rationalInterpCurve(x.map(function(y) {
+			return y.point;
+		}), x.length - 1);
+	}).map(cd => new verb.geom.NurbsCurve(cd));
+}
diff --git a/web/app/brep/operations/boolean.js b/web/app/brep/operations/boolean.js
index 4d3b5ca5..601ef77c 100644
--- a/web/app/brep/operations/boolean.js
+++ b/web/app/brep/operations/boolean.js
@@ -1,4 +1,3 @@
-import * as BREPBuilder from '../brep-builder';
 import {BREPValidator} from '../brep-validator';
 import {Edge} from '../topo/edge';
 import {Loop} from '../topo/loop';
@@ -15,10 +14,7 @@ export const TOLERANCE_HALF = TOLERANCE * 0.5;
 const DEBUG = {
   OPERANDS_MODE: false,
   LOOP_DETECTION: true,
-  FACE_FACE_INTERSECTION: false,
-  FACE_EDGE_INTERSECTION: false,
-  SEWING: false,
-  EDGE_MERGING: true,
+  FACE_FACE_INTERSECTION: true,
   NOOP: () => {}
 };
 
@@ -62,11 +58,11 @@ export function invert( shell ) {
 
 export function BooleanAlgorithm( shell1, shell2, type ) {
 
-  POINT_TO_VERT.clear();
-
   let facesData = [];
 
   mergeVertices(shell1, shell2);
+  initVertexFactory(shell1, shell2)
+
   intersectEdges(shell1, shell2);
 
   initSolveData(shell1, facesData);
@@ -75,7 +71,7 @@ export function BooleanAlgorithm( shell1, shell2, type ) {
   intersectFaces(shell1, shell2, type);
 
   for (let faceData of facesData) {
-    initGraph(faceData);
+    faceData.initGraph();
   }
 
   const allFaces = [];
@@ -119,6 +115,7 @@ function detectLoops(face) {
   const seen = new Set();
   let edges = [];
   for (let e of face.edges) edges.push(e);
+  for (let e of faceData.loopOfNew.halfEdges) edges.push(e);
   while (true) {
     let edge = edges.pop();
     if (!edge) {
@@ -127,8 +124,7 @@ function detectLoops(face) {
     if (seen.has(edge)) {
       continue;
     }
-    const loop = new Loop();
-    loop.face = face;
+    const loop = new Loop(null);
     let surface = face.surface;
     while (edge) {
       if (DEBUG.LOOP_DETECTION) {
@@ -147,24 +143,13 @@ function detectLoops(face) {
     }
 
     if (loop.halfEdges[0].vertexA === loop.halfEdges[loop.halfEdges.length - 1].vertexB) {
-      for (let halfEdge of loop.halfEdges) {
-        halfEdge.loop = loop;
-      }
-
-      BREPBuilder.linkSegments(loop.halfEdges);
+      loop.link();
       loops.push(loop);
     }
   }
   return loops;
 }
 
-function initGraph(faceData) {
-  faceData.vertexToEdge.clear();
-  for (let he of faceData.face.edges) {
-    addToListInMap(faceData.vertexToEdge, he.vertexA, he);
-  }
-}
-
 function edgeV(edge) {
   return edge.vertexB.point.minus(edge.vertexA.point)._normalize();
 }
@@ -173,7 +158,7 @@ export function mergeVertices(shell1, shell2) {
   const toSwap = new Map();
   for (let v1 of shell1.vertices) {
     for (let v2 of shell2.vertices) {
-      if (math.areVectorsEqual(v1.point, v2.point, TOLERANCE)) {
+      if (veq(v1.point, v2.point)) {
         toSwap.set(v2, v1);
       }
     }
@@ -193,13 +178,6 @@ export function mergeVertices(shell1, shell2) {
   }
 }
 
-function squash(face, edges) {
-  face.outerLoop = new Loop();
-  face.outerLoop.face = face;
-  edges.forEach(he => face.outerLoop.halfEdges.push(he));
-  face.innerLoops = [];
-}
-
 function filterFaces(faces, newLoops) {
   const validFaces = new Set(faces);
   const result = new Set();
@@ -248,6 +226,7 @@ function traverseFaces(face, validFaces, callback) {
 export function loopsToFaces(originFace, loops, out) {
   const face = new Face(originFace.surface);
   face.innerLoops = loops;
+  loops.forEach(loop => loop.face = face);
   out.push(face);
 }
 
@@ -321,16 +300,16 @@ function intersectEdges(shell1, shell2) {
           for (let point of points) {
             const {u0, u1} = point;
             let vertex;
-            if (equal(u0, 0)) {
+            if (eq(u0, 0)) {
               vertex = e1.halfEdge1.vertexA;
-            } else if (equal(u0, 1)) {
+            } else if (eq(u0, 1)) {
               vertex = e1.halfEdge1.vertexB;
-            } else if (equal(u1, 0)) {
+            } else if (eq(u1, 0)) {
               vertex = e2.halfEdge1.vertexA;
-            } else if (equal(u1, 1)) {
+            } else if (eq(u1, 1)) {
               vertex = e2.halfEdge1.vertexB;
             } else {
-              vertex = newVertex(e1.curve.point(u0));
+              vertex = vertexFactory.create(e1.curve.point(u0));
             }
             const new1 = splitEdgeByVertex(e1, vertex);
             const new2 = splitEdgeByVertex(e2, vertex);
@@ -362,7 +341,7 @@ function fixCurveDirection(curve, surface1, surface2, operationType) {
     expectedDirection._negate();
   }
   let sameAsExpected = expectedDirection.dot(tangent) > 0;
-  if (sameAsExpected) {
+  if (!sameAsExpected) {
     curve = curve.invert();
   }
   return curve;
@@ -373,7 +352,7 @@ function newEdgeDirectionValidityTest(e, curve) {
   let point = e.halfEdge1.vertexA.point;
   let tangent = curve.tangentAtPoint(point);
   assert('tangent of originated curve and first halfEdge should be the same', math.vectorsEqual(tangent, e.halfEdge1.tangent(point)));
-  assert('tangent of originated curve and second halfEdge should be the opposite', math.vectorsEqual(tangent, e.halfEdge2.tangent(point)));
+  assert('tangent of originated curve and second halfEdge should be the opposite', math.vectorsEqual(tangent._negate(), e.halfEdge2.tangent(point)));
 }
 
 function intersectFaces(shell1, shell2, operationType) {
@@ -381,13 +360,15 @@ function intersectFaces(shell1, shell2, operationType) {
   for (let i = 0; i < shell1.faces.length; i++) {
     const face1 = shell1.faces[i];
     if (DEBUG.FACE_FACE_INTERSECTION) {
-      __DEBUG__.Clear(); __DEBUG__.AddFace(face1, 0x00ff00);
+      __DEBUG__.Clear(); 
+      __DEBUG__.AddFace(face1, 0x00ff00);
       DEBUG.NOOP();
     }
     for (let j = 0; j < shell2.faces.length; j++) {
       const face2 = shell2.faces[j];
       if (DEBUG.FACE_FACE_INTERSECTION) {
-        __DEBUG__.Clear(); __DEBUG__.AddFace(face1, 0x00ff00);
+        __DEBUG__.Clear();
+        __DEBUG__.AddFace(face1, 0x00ff00);
         __DEBUG__.AddFace(face2, 0x0000ff);
         if (face1.refId === 0 && face2.refId === 0) {
           DEBUG.NOOP();
@@ -419,7 +400,7 @@ function intersectFaces(shell1, shell2, operationType) {
 
 function addNewEdge(face, halfEdge) {
   const data = face.data[MY];
-  data.newEdges.push(halfEdge);
+  data.loopOfNew.halfEdges.push(halfEdge);
   halfEdge.loop = data.loopOfNew;
   EdgeSolveData.createIfEmpty(halfEdge).newEdgeFlag = true;
   //addToListInMap(data.vertexToEdge, halfEdge.vertexA, halfEdge);
@@ -445,7 +426,7 @@ function filterNodes(nodes) {
       if (i === j) continue;
       const node2 = nodes[j];
       if (node2 !== null) {
-        if (equal(node2.u, node1.u)) {
+        if (eq(node2.u, node1.u)) {
           if (node1.normal + node2.normal === 0) {
             nodes[i] = null
           }
@@ -474,12 +455,12 @@ function intersectCurveWithEdge(curve, edge, result) {
     const {u0, u1} = point;
 
     let vertex;
-    if (equal(u0, 0)) {
+    if (eq(u0, 0)) {
       vertex = edge.edge.halfEdge1.vertexA;
-    } else if (equal(u0, 1)) {
+    } else if (eq(u0, 1)) {
       vertex = edge.edge.halfEdge1.vertexB;
     } else {
-      vertex = new Vertex(edge.edge.curve.point(u0));
+      vertex = vertexFactory.create(point.p0);
     }
 
     result.push(new Node(vertex, edge, curve, u1));
@@ -498,8 +479,8 @@ function split(nodes, curve, result) {
 
     const edge = new Edge(curve, inNode.vertex, outNode.vertex);
 
-    splitEdgeByVertex(inNode.edge, edge.halfEdge1.vertexA);
-    splitEdgeByVertex(outNode.edge, edge.halfEdge1.vertexB);
+    splitEdgeByVertex(inNode.edge.edge, edge.halfEdge1.vertexA);
+    splitEdgeByVertex(outNode.edge.edge, edge.halfEdge1.vertexB);
 
     result.push(edge);
   }
@@ -533,23 +514,15 @@ function splitEdgeByVertex(edge, vertex) {
   return [edge1, edge2];
 }
 
-const POINT_TO_VERT = new Map();
-function newVertex(point) {
-  let vertex = POINT_TO_VERT.get(point);
-  if (!vertex) {
-    vertex = new Vertex(point);
-    duplicatePointTest(point);
-    POINT_TO_VERT.set(point, vertex);
-  }
-  return vertex;
-}
-
 function nodeNormal(point, edge, curve) {
+  const normal = edge.loop.face.surface.normal(point);
   const edgeTangent = edge.tangent(point);
   const curveTangent = curve.tangentAtPoint(point);
-
-  let dot = edgeTangent.dot(curveTangent);
-  if (equal(dot, 0)) {
+  
+    
+  let cross = normal.cross(edgeTangent);
+  let dot = cross.dot(curveTangent);
+  if (eq(dot, 0)) {
     dot = 0;
   } else {
     if (dot < 0) 
@@ -597,23 +570,44 @@ function Node(vertex, edge, curve, u) {
 }
 
 
-let __DEBUG_POINT_DUPS = [];
-function duplicatePointTest(point, data) {
-  data = data || {};
-  let res = false;
-  for (let entry of __DEBUG_POINT_DUPS) {
-    let other = entry[0];
-    if (math.areVectorsEqual(point, other, TOLERANCE)) {
-      res = true;
-      break;
+let vertexFactory = null;
+function initVertexFactory(shell1, shell2) {
+  vertexFactory = new VertexFactory();
+  vertexFactory.addVertices(shell1.vertices);
+  vertexFactory.addVertices(shell2.vertices);
+}
+
+class VertexFactory {
+
+  constructor() {
+    this.vertices = [];
+  }
+
+  addVertices(vertices) {  
+    for (let v of vertices) {
+      this.vertices.push(v);
     }
   }
-  __DEBUG_POINT_DUPS.push([point, data]);
-  if (res) {
-    __DEBUG__.AddPoint(point);
-    console.error('DUPLICATE DETECTED: ' + point)
+
+  find(point) {
+    for (let vertex of this.vertices) {
+      if (veq(point, vertex.point)) {
+        return vertex;
+      }
+    }
+    return null;  
+  }
+
+  create(point) {
+    
+    let vertex = this.find(point);
+    if (vertex === null) {
+      vertex = new Vertex(point);
+      this.vertices.push(vertex);
+      console.log("DUPLICATE DETECTED: " + vertex);
+    }
+    return vertex;
   }
-  return res;
 }
 
 class SolveData {
@@ -625,11 +619,22 @@ class SolveData {
 class FaceSolveData {
   constructor(face) {
     this.face = face;
-    this.loopOfNew = new Loop();
-    this.newEdges = this.loopOfNew.halfEdges;
+    this.loopOfNew = new Loop(face);
     this.vertexToEdge = new Map();
-    this.overlaps = new Set();
-    this.loopOfNew.face = face;
+  }
+
+  initGraph() {
+    this.vertexToEdge.clear();
+    for (let he of this.face.edges) {
+      this.addToGraph(he);
+    }
+    for (let he of this.loopOfNew.halfEdges) {
+      this.addToGraph(he);
+    }
+  }
+
+  addToGraph(he) {
+    addToListInMap(this.vertexToEdge, he.vertexA, he);
   }
 }
 
@@ -660,10 +665,14 @@ function __DEBUG_OPERANDS(shell1, shell2) {
   }
 }
 
-function equal(v1, v2) {
+function eq(v1, v2) {
   return math.areEqual(v1, v2, TOLERANCE);
 }
 
+function veq(v1, v2) {
+  return math.areVectorsEqual(v1, v2, TOLERANCE);
+}
+
 function assert(name, cond) {
   if (!cond) {
     throw 'ASSERTION FAILED: ' + name;
diff --git a/web/app/brep/topo/face.js b/web/app/brep/topo/face.js
index 5c9602a1..64aab3d1 100644
--- a/web/app/brep/topo/face.js
+++ b/web/app/brep/topo/face.js
@@ -8,7 +8,7 @@ export class Face extends TopoObject {
     this.id = undefined;
     this.surface = surface;
     this.shell = null;
-    this.outerLoop = new Loop();
+    this.outerLoop = new Loop(this);
     this.innerLoops = [];
     this.defineIterable('loops', () => loopsGenerator(this));
     this.defineIterable('edges', () => halfEdgesGenerator(this))
diff --git a/web/app/brep/topo/loop.js b/web/app/brep/topo/loop.js
index 74b62623..c28dd7b3 100644
--- a/web/app/brep/topo/loop.js
+++ b/web/app/brep/topo/loop.js
@@ -4,9 +4,9 @@ import * as math from '../../math/math'
 
 export class Loop extends TopoObject {
 
-  constructor() {
+  constructor(face) {
     super();
-    this.face = null;
+    this.face = face;
     this.halfEdges = [];
   }
 
@@ -26,6 +26,8 @@ export class Loop extends TopoObject {
       const next = this.halfEdges[j];
       curr.next = next;
       next.prev = curr;
+
+      curr.loop = this;
     }
   }
 }
diff --git a/web/app/math/math.js b/web/app/math/math.js
index 6dc0c5b3..15abe4ad 100644
--- a/web/app/math/math.js
+++ b/web/app/math/math.js
@@ -58,6 +58,12 @@ export function areVectorsEqual(v1, v2, tolerance) {
          areEqual(v1.z, v2.z, tolerance);
 }
 
+export function areVectorsEqual3(v1, v2, tolerance) {
+  return areEqual(v1[0], v2[0], tolerance) &&
+         areEqual(v1[1], v2[1], tolerance) &&
+         areEqual(v1[2], v2[2], tolerance);
+}
+
 export function vectorsEqual(v1, v2) {
   return areVectorsEqual(v1, v2, TOLERANCE);
 }