diff --git a/web/app/3d/modeler-app.js b/web/app/3d/modeler-app.js
index 7e10bff0..abad23b9 100644
--- a/web/app/3d/modeler-app.js
+++ b/web/app/3d/modeler-app.js
@@ -92,27 +92,32 @@ App.prototype.scratchCode = function() {
 
   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);
+  let result = app.TPI.brep.bool.subtract(box1, box2);
+  result = app.TPI.brep.bool.subtract(result, box3);
   // app.addShellOnScene(box1);
-  // app.addShellOnScene(result);
+  app.addShellOnScene(result);
 
 
-  let curve1 = new NurbsCurve(new verb.geom.NurbsCurve({"degree":6,"controlPoints":[[150,149.99999999999997,-249.99999999999994,1],[108.33333333333051,150.00000000000907,-250.00000000001975,1],[66.6666666666712,149.99999999998562,-249.99999999996987,1],[24.99999999999545,150.00000000001364,-250.00000000002711,1],[-16.66666666666362,149.99999999999145,-249.9999999999837,1],[-58.33333333333436,150.0000000000029,-250.00000000000531,1],[-99.99999999999997,150,-250,1]],"knots":[0,0,0,0,0,0,0,1,1,1,1,1,1,1]}));
-  let curve2 = new NurbsCurve(new verb.geom.NurbsCurve({"degree":9,"controlPoints":[[100,-250,-250,1],[99.9999999999927,-194.44444444444687,-250.00000000000028,1],[100.00000000002228,-138.8888888888811,-249.99999999999838,1],[99.99999999995923,-83.33333333334777,-250.00000000000287,1],[100.00000000005268,-27.77777777775936,-249.99999999999744,1],[99.9999999999493,27.777777777760704,-250.0000000000008,1],[100.00000000003591,83.33333333334477,-250.00000000000063,1],[99.99999999998269,138.88888888888374,-249.99999999999966,1],[100.00000000000443,194.44444444444562,-249.99999999999986,1],[100,250,-250,1]],"knots":[0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1]}));
+  // let curve1 = new NurbsCurve(new verb.geom.NurbsCurve({"degree":6,"controlPoints":[[150,149.99999999999997,-249.99999999999994,1],[108.33333333333051,150.00000000000907,-250.00000000001975,1],[66.6666666666712,149.99999999998562,-249.99999999996987,1],[24.99999999999545,150.00000000001364,-250.00000000002711,1],[-16.66666666666362,149.99999999999145,-249.9999999999837,1],[-58.33333333333436,150.0000000000029,-250.00000000000531,1],[-99.99999999999997,150,-250,1]],"knots":[0,0,0,0,0,0,0,1,1,1,1,1,1,1]}));
+  // let curve2 = new NurbsCurve(new verb.geom.NurbsCurve({"degree":9,"controlPoints":[[100,-250,-250,1],[99.9999999999927,-194.44444444444687,-250.00000000000028,1],[100.00000000002228,-138.8888888888811,-249.99999999999838,1],[99.99999999995923,-83.33333333334777,-250.00000000000287,1],[100.00000000005268,-27.77777777775936,-249.99999999999744,1],[99.9999999999493,27.777777777760704,-250.0000000000008,1],[100.00000000003591,83.33333333334477,-250.00000000000063,1],[99.99999999998269,138.88888888888374,-249.99999999999966,1],[100.00000000000443,194.44444444444562,-249.99999999999986,1],[100,250,-250,1]],"knots":[0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1]}));
 
-  __DEBUG__.AddCurve(curve1);
-  __DEBUG__.AddCurve(curve2);
+  // var p1 = [-50,0,0], p2 = [100,0,0], p3 = [100,100,0], p4 = [0,100,0], p5 = [50, 50, 0];
+  // var pts = [p1, p2, p3, p4, p5];
+  // let curve1 = new NurbsCurve(verb.geom.NurbsCurve.byPoints( pts, 3 ));
+  //
+  // var p1a = [-50,0,0], p2a = [50,-10,0], p3a = [150,50,0], p4a = [30,100,0], p5a = [50, 120, 0];
+  // var ptsa = [p1a, p2a, p3a, p4a, p5a];
+  // let curve2 = new NurbsCurve(verb.geom.NurbsCurve.byPoints( ptsa, 3 ));
+  //
+  // __DEBUG__.AddCurve(curve1);
+  // __DEBUG__.AddCurve(curve2);
 
 
-  // let curves =  surface.intersectSurface(box1.faces[0].surface);
-  // const curve = box1.faces[0].outerLoop.halfEdges[0].edge.curve;
-  let points = curve1.intersectCurve(curve2);
-  for (let p of points) {
-    __DEBUG__.AddPoint(p.p0);
-  }
-
 
+  // let points = curve1.intersectCurve(curve2);
+  // for (let p of points) {
+  //   __DEBUG__.AddPoint(p.p0);
+  // }
 
   app.viewer.render();
 };
diff --git a/web/app/brep/geom/impl/nurbs.js b/web/app/brep/geom/impl/nurbs.js
index da86def0..42bce564 100644
--- a/web/app/brep/geom/impl/nurbs.js
+++ b/web/app/brep/geom/impl/nurbs.js
@@ -3,6 +3,7 @@ import * as math from  '../../../math/math'
 import {Point} from '../point'
 import {Surface} from "../surface";
 import Vector from "../../../math/vector";
+import * as vec from "../../../math/vec";
 import {Curve} from "../curve";
 
 export class NurbsCurve extends Curve {
@@ -76,7 +77,7 @@ export class NurbsCurve extends Curve {
     isecOn(other, this, 0);
     isecOn(other, this, 1);
 
-    verb.geom.Intersect.curves(this.verb, other.verb, tol).forEach( i => add({
+    verb_curve_isec(this.verb, other.verb, tol).forEach( i => add({
       u0: i.u0,
       u1: i.u1,
       p0: i.point0,
@@ -86,10 +87,11 @@ export class NurbsCurve extends Curve {
       i.p0 = pt(i.p0);
       i.p1 = pt(i.p1);
     });
-    return isecs.filter(({u0, u1}) => {
-      return Math.abs(this.tangentAtParam(u0).dot(other.tangentAtParam(u1))) <= tol;
+    isecs = isecs.filter(({u0, u1}) => {
+      let collinearFactor = Math.abs(this.tangentAtParam(u0).dot(other.tangentAtParam(u1)));
+      return !math.areEqual(collinearFactor, 1, tol);
     });
-
+    return isecs;
 }
 
   static createByPoints(points, degeree) {
@@ -146,7 +148,7 @@ export class NurbsSurface extends Surface {
   }
 
   intersectSurfaceForSameClass(other, tol) {
-    const curves = verb_isec(this.verb, other.verb);
+    const curves = verb_surface_isec(this.verb, other.verb, tol);
     let inverted = this.inverted !== other.inverted;
     return curves.map(curve => new NurbsCurve(inverted ?  curve.reverse() : curve));
   }
@@ -180,15 +182,13 @@ function pt(data) {
   return new Point().set3(data);
 }
 
-
-function verb_isec(nurbs1, nurbs2) {
-  const tol = 1e-3
+function verb_surface_isec(nurbs1, nurbs2, tol) {
   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) {
+	const tess1 = verb.eval.Tess.rationalSurfaceAdaptive(surface0);
+	const tess2 = verb.eval.Tess.rationalSurfaceAdaptive(surface1);
+	const resApprox = verb.eval.Intersect.meshes(tess1,tess2);
+	const exactPls = resApprox.map(function(pl) {
 		return pl.map(function(inter) {
 			return verb.eval.Intersect.surfacesAtPointWithEstimate(surface0,surface1,inter.uv0,inter.uv1,tol);
 		});
@@ -199,3 +199,70 @@ function verb_isec(nurbs1, nurbs2) {
 		}), x.length - 1);
 	}).map(cd => new verb.geom.NurbsCurve(cd));
 }
+
+function verb_curve_isec(curve1, curve2, tol) {
+
+  let result = [];
+  let segs1 = curve1.tessellate();
+  let segs2 = curve2.tessellate();
+
+  for (let i = 0; i < segs1.length - 1; i++) {
+    let a1 = segs1[i];
+    let b1 = segs1[i + 1];
+    for (let j = 0; j < segs2.length - 1; j++) {
+      let a2 = segs2[j];
+      let b2 = segs2[j + 1];
+
+      //TODO: minimize
+      let isec = intersectSegs(a1, b1, a2, b2, tol);
+      if (isec !== null) {
+        let {u1, u2, point1, point2, l1, l2} = isec;
+        result.push({
+          u0: curve1.closestParam(point1),
+          u1: curve2.closestParam(point2),
+          point0: point1,
+          point1: point2
+        });
+        if (math.areEqual(u1, l1, tol )) {
+          i ++;
+        }
+        if (math.areEqual(u2, l2, tol )) {
+          j ++;
+        }
+      }
+    }
+  }
+  return result;
+}
+
+export function lineLineIntersection(p1, p2, v1, v2) {
+  let zAx = vec.cross(v1, v2);
+  const n1 = vec._normalize(vec.cross(zAx, v1));
+  const n2 = vec._normalize(vec.cross(zAx, v2));
+  return {
+    u1: vec.dot(n2, vec.sub(p2, p1)) / vec.dot(n2, v1),
+    u2: vec.dot(n1, vec.sub(p1, p2)) / vec.dot(n1, v2),
+  }
+}
+
+function intersectSegs(a1, b1, a2, b2, tol) {
+  let v1 = vec.sub(b1, a1);
+  let v2 = vec.sub(b2, a2);
+  let l1 = vec.length(v1);
+  let l2 = vec.length(v2);
+  vec._div(v1, l1);
+  vec._div(v2, l2);
+
+  let {u1, u2} = lineLineIntersection(a1, a2, v1, v2);
+  let point1 = vec.add(a1, vec.mul(v1, u1));
+  let point2 = vec.add(a2, vec.mul(v2, u2));
+  let p2p = vec.lengthSq(vec.sub(point1, point2));
+  let eq = (a, b) => math.areEqual(a, b, tol);
+  if (u1 !== Infinity && u2 !== Infinity && math.areEqual(p2p, 0, tol*tol) &&
+      ((u1 >0 && u1 < l1) || eq(u1, 0) || eq(u1, l1)) &&
+      ((u2 >0 && u2 < l2) || eq(u2, 0) || eq(u2, l2))
+     ) {
+    return {point1, point2, u1, u2, l1, l2}
+  }
+  return null;
+}
diff --git a/web/app/brep/operations/boolean.js b/web/app/brep/operations/boolean.js
index 124f99c2..9463f0ed 100644
--- a/web/app/brep/operations/boolean.js
+++ b/web/app/brep/operations/boolean.js
@@ -7,8 +7,8 @@ import {Vertex} from '../topo/vertex';
 import Vector from '../../math/vector';
 import * as math from '../../math/math';
 
-export const TOLERANCE = 1e-8;
- export const TOLERANCE_SQ = TOLERANCE * TOLERANCE;
+export const TOLERANCE = 1e-3;
+export const TOLERANCE_SQ = TOLERANCE * TOLERANCE;
 export const TOLERANCE_HALF = TOLERANCE * 0.5;
 
 const DEBUG = {
@@ -448,6 +448,7 @@ function collectNodesOfIntersection(curve, loop, nodes) {
 function intersectCurveWithEdge(curve, edge, result) {
   __DEBUG__.AddCurve(curve, 0xffffff);
   __DEBUG__.AddHalfEdge(edge, 0xff00ff);
+  console.log(1)
   const points = edge.edge.curve.intersectCurve(curve, TOLERANCE);
   for (let point of points) {
     const {u0, u1} = point;
diff --git a/web/app/math/vec.js b/web/app/math/vec.js
new file mode 100644
index 00000000..0ebc9feb
--- /dev/null
+++ b/web/app/math/vec.js
@@ -0,0 +1,122 @@
+export function __mul(v, scalar, out) {
+  out[0] = scalar * v[0];
+  out[1] = scalar * v[1];
+  out[2] = scalar * v[2];
+  return out;
+}
+
+export function _mul(v, scalar) {
+  return __mul(v, scalar, v);
+}
+
+export function mul(v, scalar) {
+  return __mul(v, scalar, []);
+}
+
+export function __div(v, scalar, out) {
+  out[0] = v[0] / scalar;
+  out[1] = v[1] / scalar;
+  out[2] = v[2] / scalar;
+  return out;
+}
+
+export function _div(v, scalar) {
+  return __div(v, scalar, v);
+}
+
+export function div(v, scalar) {
+  return __div(v, scalar, []);
+}
+
+
+export function __add(v1, v2, out) {
+  out[0] = v1[0] + v2[0];
+  out[1] = v1[1] + v2[1];
+  out[2] = v1[2] + v2[2];
+  return out;
+}
+
+export function _add(v1, v2) {
+  return __add(v1, v2, v1);
+}
+
+export function add(v1, v2) {
+  return __add(v1, v2, []);
+}
+
+export function __sub(v1, v2, out) {
+  out[0] = v1[0] - v2[0];
+  out[1] = v1[1] - v2[1];
+  out[2] = v1[2] - v2[2];
+  return out;
+}
+
+export function _sub(v1, v2) {
+  return __sub(v1, v2, v1);
+}
+
+export function sub(v1, v2) {
+  return __sub(v1, v2, []);
+}
+
+
+export function __negate(v, out) {
+  out[0] = - v[0];
+  out[1] = - v[1];
+  out[2] = - v[2];
+  return out;
+}
+
+export function _negate(v) {
+  return __negate(v, v);
+}
+
+export function negate(v) {
+  return __negate(v, []);
+}
+
+
+export function dot(v1, v2) {
+  return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
+}
+
+export function __cross(v1, v2, out) {
+  out[0] = v1[1] * v2[2] - v1[2] * v2[1];
+  out[1] = v1[2] * v2[0] - v1[0] * v2[2];
+  out[2] = v1[0] * v2[1] - v1[1] * v2[0];
+  return out;
+}
+
+export function _cross(v1, v2) {
+  return __cross(v1, v2);
+}
+
+export function cross(v1, v2) {
+  return __cross(v1, v2, []);
+}
+
+
+export function __normalize(v, out) {
+  const mag = length(v);
+  if (mag === 0.0) {
+    out[0] = out[1] = out[2] = 0;
+  }
+  return __div(v, mag, out)
+}
+
+export function _normalize(v) {
+  return __normalize(v, v);
+}
+
+export function normalize(v) {
+  return __normalize(v, []);
+}
+
+
+export function lengthSq(v) {
+  return dot(v, v);
+}
+
+export function length(v) {
+  return Math.sqrt(lengthSq(v));
+}