diff --git a/package.json b/package.json
index 0ed9e0ab..7d91f714 100644
--- a/package.json
+++ b/package.json
@@ -42,6 +42,7 @@
     "grunt-contrib-copy": "1.0.0"
   },
   "dependencies": {
+    "clipper-lib": "6.2.1",
     "diff-match-patch": "1.0.0",
     "earcut": "2.1.1",
     "handlebars": "4.0.6",
diff --git a/web/app/3d/debug.js b/web/app/3d/debug.js
index f30492b0..86166bb0 100644
--- a/web/app/3d/debug.js
+++ b/web/app/3d/debug.js
@@ -56,6 +56,15 @@ function addGlobalDebugActions(app) {
       }
     },
 
+    AddPointPolygons: (polygons, color) => {
+      for (let points of polygons) {
+        for (let i = 0; i < points.length; i ++) {
+          debugGroup.add(createLine(points[i], points[(i + 1) % points.length], color));
+        }
+      }
+      app.viewer.render();
+    },
+
     AddPlane: (plane) => {
       const geo = new THREE.PlaneBufferGeometry(2000, 2000, 8, 8);
       const coplanarPoint = plane.normal.multiply(plane.w);
diff --git a/web/app/3d/modeler-app.js b/web/app/3d/modeler-app.js
index 03a504fc..72ad6e08 100644
--- a/web/app/3d/modeler-app.js
+++ b/web/app/3d/modeler-app.js
@@ -185,28 +185,28 @@ App.prototype.test5 = function() {
     return bb.vertex(pt.x, pt.y, pt.z);
   }
 
-  const a = vx(0.1, 0.1);
-  const b = vx(0.9, 0.1);
+  const a = vx(0.13, 0.13);
+  const b = vx(0.9, 0.13);
   const c = vx(0.9, 0.9);
-  const d = vx(0.1, 0.9);
+  const d = vx(0.13, 0.9);
 
-  const e = vx(0.3, 0.3);
-  const f = vx(0.3, 0.7);
-  const g = vx(0.7, 0.7);
-  const h = vx(0.7, 0.3);
+  const e = vx(0.33, 0.33);
+  const f = vx(0.33, 0.73);
+  const g = vx(0.73, 0.73);
+  const h = vx(0.73, 0.33);
 
 
   let shell = bb.face(srf)
     .loop()
-    .edgeTrim(a, b, srf.verb.isocurve(0.1, true))
+    .edgeTrim(a, b, srf.verb.isocurve(0.13, true))
     .edgeTrim(b, c, srf.verb.isocurve(0.9, false))
     .edgeTrim(c, d, srf.verb.isocurve(0.9, true).reverse())
-    .edgeTrim(d, a, srf.verb.isocurve(0.1, false).reverse())
+    .edgeTrim(d, a, srf.verb.isocurve(0.13, false).reverse())
     .loop()
-    .edgeTrim(e, f, srf.verb.isocurve(0.3, false))
-    .edgeTrim(f, g, srf.verb.isocurve(0.7, true))
-    .edgeTrim(g, h, srf.verb.isocurve(0.7, false).reverse())
-    .edgeTrim(h, e, srf.verb.isocurve(0.3, true).reverse())
+    .edgeTrim(e, f, srf.verb.isocurve(0.33, false))
+    .edgeTrim(f, g, srf.verb.isocurve(0.73, true))
+    .edgeTrim(g, h, srf.verb.isocurve(0.73, false).reverse())
+    .edgeTrim(h, e, srf.verb.isocurve(0.33, true).reverse())
     .build();
 
   this.addShellOnScene(shell);
@@ -214,7 +214,8 @@ App.prototype.test5 = function() {
 
 App.prototype.scratchCode = function() {
   const app = this;
-  this.cylTest();
+  this.test5();
+  // this.cylTest();
 
 
 
diff --git a/web/app/3d/tess/brep-tess.js b/web/app/3d/tess/brep-tess.js
index 2d7c33e4..51dfb824 100644
--- a/web/app/3d/tess/brep-tess.js
+++ b/web/app/3d/tess/brep-tess.js
@@ -1,6 +1,9 @@
 import PIP from "./pip";
 import earcut from 'earcut'
 import Vector from "../../math/vector";
+import {NurbsSurface} from "../../brep/geom/impl/nurbs";
+import ClipperLib from 'clipper-lib';
+import libtess from 'libtess'
 
 export default function A(face) {
 
@@ -21,122 +24,102 @@ export default function A(face) {
     }
   }
 
-  let steinerPoints = [];
   let tess = face.surface.verb.tessellate({maxDepth: 3});
-  for (let i = 0; i < tess.points.length; i++) {
-    steinerPoints.push(face.surface.createWorkingPoint(tess.uvs[i], Vector.fromData(tess.points[i])));
-  }
+  let nurbsTriangles = tess.faces.map(f => f.map(i => face.surface.createWorkingPoint(tess.uvs[i], Vector.fromData(tess.points[i]))));
 
-  let [outer, ...inners] = loops;
-  inners.forEach(inner => inner.reverse());
-  let pip = PIP(outer, inners);
-  steinerPoints = steinerPoints.filter(pt => pip(pt).inside);
+  let paths = clip(nurbsTriangles, loops);
 
-  let points = [];
-  let pointsData = [];
-  let holes = [];
+  let triangles = tessPaths(paths);
 
-  function pushLoop(loop) {
-    for (let pt of loop) {
-      pointsData.push(pt.x);
-      pointsData.push(pt.y);
-      points.push(pt);
-    }
-  }
+  let out = convertPoints(triangles, p => face.surface.workingPointTo3D(p) );
+  // __DEBUG__.AddPointPolygons(out, 0x00ffff);
+  return out;
+}
 
-  pushLoop(outer);
+function convertPoints(paths, converter) {
+  return paths.map( path => path.map(p => converter(p) ))
+}
 
-  for (let inner of inners) {
-    holes.push(pointsData.length / 2);
-    pushLoop(inner);
-  }
+function clip(triangles, loops) {
 
-  let trs = earcut(pointsData, holes);
+  const scale = 1e3 ;// multiplying by NurbsSurface.WORKING_POINT_SCALE_FACTOR gives 1e6
 
-  let triangles = [];
 
-  for (let i = 0; i < trs.length; i += 3) {
-    const tr = [trs[i], trs[i + 1], trs[i + 2]];
+  let clip_paths = convertPoints(loops, p => ({X:p.x, Y:p.y}) );
+  ClipperLib.JS.ScaleUpPaths(clip_paths, scale);
 
-    __DEBUG__.AddPointPolygon(tr.map( ii => new Vector(pointsData[ii * 2], pointsData[ii * 2 + 1], 0) ));
-
-    triangles.push(tr.map(i => points[i]));
-  }
-
-  splitTriangles(triangles, steinerPoints);
-
-  triangles = triangles.filter(tr => tr !== null);
+  let out = [];
 
   for (let tr of triangles) {
-    for (let i = 0; i < tr.length; i++) {
-      tr[i] = tr[i].__3D;
-    }
+    let cpr = new ClipperLib.Clipper();
+    let subj_paths  = convertPoints([tr], p => ({X:p.x, Y:p.y}) );
+
+    ClipperLib.JS.ScaleUpPaths(subj_paths, scale);
+
+
+    cpr.AddPaths(subj_paths, ClipperLib.PolyType.ptSubject, true);  // true means closed path
+    cpr.AddPaths(clip_paths, ClipperLib.PolyType.ptClip, true);
+
+    let solution_paths = new ClipperLib.Paths();
+    let succeeded = cpr.Execute(ClipperLib.ClipType.ctIntersection, solution_paths, ClipperLib.PolyFillType.pftNonZero, ClipperLib.PolyFillType.pftNonZero);
+    ClipperLib.JS.ScaleUpPaths(solution_paths, 1.0/scale);
+    solution_paths.forEach(p => out.push(p));
   }
 
-  return triangles;
+  out = convertPoints(out, p => new Vector(p.X, p.Y, 0) );
+  return out;
 }
 
-function splitTriangles(triangles, steinerPoints) {
-  for (let sp of steinerPoints) {
-    __DEBUG__.AddPoint(sp);
-    let newTrs = [];
-    for (let i = 0; i < triangles.length; ++i) {
-      let tr = triangles[i];
-      if (tr === null) {
-        continue;
-      }
-      let pip = new PIP(tr);
-      let res = pip(sp);
-      if (!res.inside || res.vertex) {
-        continue;
-      } else {
-        if (res.edge) {
-          let [tr1, tr2] = splitEdgeOfTriangle(sp, tr, res.edge);
-          if (tr1 && tr2) {
-            newTrs.push(tr1, tr2);
-            triangles[i] = null;
-          }
-        } else {
-          let [tr1, tr2, tr3] = splitTriangle(sp, tr, res.edge, triangles);
-          newTrs.push(tr1, tr2, tr3);
-          triangles[i] = null;
-        }
 
-      }
+function tessPaths(paths) {
+
+  function vertexCallback(data, out) {
+    out.push(data);
+  }
+  function combinecallback(coords, data, weight) {
+  }
+  function edgeCallback(flag) {
+  }
+
+  const tessy = new libtess.GluTesselator();
+  // tessy.gluTessProperty(libtess.gluEnum.GLU_TESS_WINDING_RULE, libtess.windingRule.GLU_TESS_WINDING_POSITIVE);
+  tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_VERTEX_DATA, vertexCallback);
+  tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_BEGIN, begincallback);
+  tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_ERROR, errorcallback);
+  tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_COMBINE, combinecallback);
+  tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_EDGE_FLAG, edgeCallback);
+
+  tessy.gluTessNormal(0, 0, 1);
+
+  const vertices = [];
+  tessy.gluTessBeginPolygon(vertices);
+
+  for (let path of paths) {
+    tessy.gluTessBeginContour();
+    for (let p of path) {
+      tessy.gluTessVertex([p.x, p.y, 0], p);
     }
-    newTrs.forEach(tr => triangles.push(tr));
+    tessy.gluTessEndContour();
+  }
+  tessy.gluTessEndPolygon();
+
+  const triangled = [];
+
+  for (let i = 0;  i < vertices.length; i += 3 ) {
+    const a = vertices[i];
+    const b = vertices[i + 1];
+    const c = vertices[i + 2];
+    triangled.push([a, b, c]);
+  }
+  return triangled;
+}
+
+function begincallback(type) {
+  if (type !== libtess.primitiveType.GL_TRIANGLES) {
+    console.log('expected TRIANGLES but got type: ' + type);
   }
 }
-
-function splitEdgeOfTriangle(p, tr, edge) {
-  let n = tr.length;
-  for (let i1 = 0; i1 < n; i1 ++ ) {
-    let i2 = (i1 + 1) % n;
-    let i3 = (i1 + 2) % n;
-    if (tr[i1] === edge[0] && tr[i2] === edge[1]) {
-      let tr1 = [tr[i1], p, tr[i3]];
-      let tr2 = [p, tr[i2], tr[i3]];
-      return [tr1, tr2];
-    }
-  }
-  return [];
+function errorcallback(errno) {
+  console.log('tesselation error');
+  console.log('error number: ' + errno);
 }
-
-function splitTriangle(p, tr, edge) {
-  let [a, b, c] = tr;
-  return [
-    [a, b, p],
-    [b, c, p],
-    [c, a, p]
-  ];
-}
-
-export function isMirrored(surface) {
-  let a = surface.point(0, 0);
-  let b = surface.point(1, 0);
-  let c = surface.point(1, 1);
-  return b.minus(a).cross(c.minus(a))._normalize().dot(surface.normalUV(0, 0)) < 0;
-}
-
-
-
diff --git a/web/app/brep/brep-primitives.js b/web/app/brep/brep-primitives.js
index e9db8e3f..1152a843 100644
--- a/web/app/brep/brep-primitives.js
+++ b/web/app/brep/brep-primitives.js
@@ -21,9 +21,9 @@ export function box(w, h, d, tr) {
 
 
 export function cylinder(r, h, tr) {
-  let circle1 = new Circle(-1, new Point(0,0,0), r).toNurbs(Plane.XY);
-  let circle2 = circle1.translate(new Point(0,h,0))
-  return enclose(circle1, circle2);
+  let circle1 = new Circle(-1, new Point(0,0,h), r).toNurbs(Plane.XY);
+  let circle2 = circle1.translate(new Point(0,0,-h));
+  return enclose([circle1], [circle2]);
 }
 
 
diff --git a/web/app/brep/geom/impl/nurbs.js b/web/app/brep/geom/impl/nurbs.js
index d7ad546d..c50da49f 100644
--- a/web/app/brep/geom/impl/nurbs.js
+++ b/web/app/brep/geom/impl/nurbs.js
@@ -158,7 +158,7 @@ export class NurbsSurface extends Surface {
   }
 
   createWorkingPoint(uv, pt3d) {
-    const wp = new Vector(uv[0], uv[1], 0)._multiply(1000);
+    const wp = new Vector(uv[0], uv[1], 0)._multiply(NurbsSurface.WORKING_POINT_SCALE_FACTOR);
     if (this.mirrored) {
       wp.x *= -1;
     }
@@ -166,6 +166,17 @@ export class NurbsSurface extends Surface {
     return wp; 
   }
 
+  workingPointTo3D(wp) {
+    if (wp.__3D === undefined) {
+      const uv = wp.multiply(NurbsSurface.WORKING_POINT_UNSCALE_FACTOR);
+      if (this.mirrored) {
+        uv.x *= -1;
+      }
+      wp.__3D = this.point(uv.x, uv.y);
+    }
+    return wp.__3D;
+  }
+
   static isMirrored(surface) {  
     let a = surface.point(0, 0);
     let b = surface.point(1, 0);
@@ -198,6 +209,9 @@ export class NurbsSurface extends Surface {
   }
 }
 
+NurbsSurface.WORKING_POINT_SCALE_FACTOR = 1000;
+NurbsSurface.WORKING_POINT_UNSCALE_FACTOR = 1 / NurbsSurface.WORKING_POINT_SCALE_FACTOR;
+
 function dist(p1, p2) {
   return math.distance3(p1[0], p1[1], p1[2], p2[0], p2[1], p2[2]);
 }