diff --git a/.gitignore b/.gitignore
index 04bfedb..883286e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -214,3 +214,4 @@ __marimo__/
 
 # 
 .pyc
+__pycache__
diff --git a/examples/basic_example.ipynb b/examples/basic_example.ipynb
index a74468c..2385cd7 100644
--- a/examples/basic_example.ipynb
+++ b/examples/basic_example.ipynb
@@ -27,7 +27,7 @@
    "outputs": [],
    "source": [
     "\n",
-    "domain = box(0, 0, 200, 150)\n",
+    "domain = box(0, 0, 200, 200)\n",
     "poly1 = box(50, 50, 75, 75)\n",
     "# generate sine wave curve for river\n",
     "river_coords = [(x, 75 + 20 * np.sin(0.1 * x)) for x in range(-10, 220, 10)]\n",
@@ -36,60 +36,139 @@
     "fault_line = LineString([(100, 0), (100, 150)])"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c970674e",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "d52a10b8",
    "metadata": {},
    "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "71f205cc",
+   "metadata": {},
+   "outputs": [],
    "source": [
     "# 1. Setup Blueprint\n",
+    "background_lc=100\n",
+    "\n",
     "blueprint = ConceptualMesh(crs=\"EPSG:3857\")\n",
-    "blueprint.add_polygon(domain, zone_id=1, resolution=200)\n",
-    "blueprint.add_polygon(poly1, zone_id=2, resolution=5, z_order=1) \n",
-    "blueprint.add_point(well_point, point_id=\"Well-A\", resolution=100)\n",
-    "blueprint.add_line(fault_line, line_id=\"Fault-1\", resolution=5, \n",
-    "                   is_barrier=True)\n",
-    "blueprint.add_line(river_coords, line_id=\"riv-1\", resolution=3, \n",
-    "                   is_barrier=False)\n",
+    "blueprint.add_polygon(domain, zone_id=1)#,border_density=100,dist_max_in=background_lc)\n",
+    "blueprint.add_polygon(poly1, zone_id=2, \n",
+    "                      resolution=1, z_order=2,\n",
+    "                      dist_max_out=3* background_lc,\n",
+    "                      #dist_max_in=10\n",
+    "                      ) \n",
+    "blueprint.add_point(well_point, point_id=\"Well-A\", \n",
+    "                    resolution=2,\n",
+    "                    dist_max=300)\n",
+    "blueprint.add_line(fault_line, line_id=\"Fault-1\", \n",
+    "                   resolution=1, \n",
+    "                   is_barrier=True,\n",
+    "                   #straddle_width=1,\n",
+    "                   dist_max=3* 5)\n",
+    "blueprint.add_line(river_coords, line_id=\"riv-1\", \n",
+    "                   resolution=1, \n",
+    "                  is_barrier=False,\n",
+    "                  dist_max=3* 5,)\n",
     "\n",
     "clean_polys, clean_lines, clean_pts = blueprint.generate()\n",
     "\n",
     "# 2. Mesh Generation\n",
-    "mesher = MeshGenerator()\n",
+    "mesher = MeshGenerator(background_lc=background_lc, verbosity=10)\n",
     "mesher.generate(clean_polys, clean_lines, clean_pts)\n",
     "\n",
     "# 3. Voronoi Conversion\n",
-    "tessellator = VoronoiTessellator(mesher, blueprint)\n",
+    "tessellator = VoronoiTessellator(mesher, blueprint, clip_to_boundary=True)\n",
     "grid_gdf = tessellator.generate()\n",
     "\n",
     "# 4. Output\n",
-    "#grid_gdf.to_file(\"mf6_grid.shp\")"
+    "#grid_gdf.to_file(\"mf6_grid.shp\")\n",
+    "\n",
+    "\n",
+    "\n",
+    "fig,ax = plt.subplots(1,1,figsize=(12,12))\n",
+    "ax.set_aspect('equal')\n",
+    "# plot domain\n",
+    "ax.plot(*domain.exterior.xy,marker='.', color='black')\n",
+    "# plot poly1\n",
+    "#ax.plot(*poly1.exterior.xy, color='blue',lw=1,ls='--',label='Polygon 1')\n",
+    "ax.plot(*river_coords.xy, color='cyan',lw=1,ls='--',label='River')\n",
+    "ax.plot(*fault_line.xy, color='r',lw=1,ls='--',label='Fault')\n",
+    "ax.scatter(*well_point.xy, color='r',s=20,marker='x',label='Well')\n",
+    "grid_gdf.plot(ax=ax, alpha=0.5, edgecolor='k', linewidth=0.2)\n",
+    "\n",
+    "ax.legend(loc='upper left', fontsize=8)\n",
+    "#ax.set_ylim(50,80)\n",
+    "#ax.set_xlim(90, 110)\n",
+    "#\n",
+    "#ax.set_ylim(40,80)\n",
+    "#ax.set_xlim(40, 80)\n",
+    "\n",
+    "fig.tight_layout()"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "71f205cc",
+   "id": "efd5f0a2",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "eda65ca1",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0ba0153a",
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig,ax = plt.subplots(1,1,figsize=(15,15))\n",
-    "ax.set_aspect('equal')\n",
-    "# plot domain\n",
-    "ax.plot(*domain.exterior.xy,marker='.', color='black')\n",
-    "# plot poly1\n",
-    "ax.plot(*poly1.exterior.xy, color='blue',lw=1,ls='--')\n",
-    "ax.plot(*river_coords.xy, color='cyan',lw=1,ls='--')\n",
-    "ax.plot(*fault_line.xy, color='r',lw=1,ls='--')\n",
-    "ax.scatter(*well_point.xy, color='green',s=100,marker='x')\n",
-    "grid_gdf.plot(ax=ax)"
+    "from vorflow.utils import *"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cfd21be8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gdf = calculate_mesh_quality(grid_gdf,calc_ortho=True)\n",
+    "gdf.plot(column='drift_ratio', cmap='viridis', legend=True, figsize=(10,8),vmax=.25)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9f7361dd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "summarize_quality(gdf)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "b6c2a945",
+   "id": "870ec1be",
    "metadata": {},
    "outputs": [],
    "source": []
diff --git a/examples/comprehensive_demo.ipynb b/examples/comprehensive_demo.ipynb
new file mode 100644
index 0000000..8df3770
--- /dev/null
+++ b/examples/comprehensive_demo.ipynb
@@ -0,0 +1,266 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "6f3f8101",
+   "metadata": {},
+   "source": [
+    "# Vorflow Comprehensive Demo\n",
+    "\n",
+    "This notebook demonstrates the full capabilities of the `vorflow` package for generating unstructured Voronoi grids for MODFLOW 6.\n",
+    "\n",
+    "## Workflow\n",
+    "1. **Define Geometry**: Create Shapely objects (Polygons, Lines, Points).\n",
+    "2. **ConceptualMesh**: Register these objects with refinement parameters.\n",
+    "3. **MeshGenerator**: Generate a triangular mesh using Gmsh.\n",
+    "4. **VoronoiTessellator**: Convert the triangles to a Voronoi grid.\n",
+    "5. **Analysis**: Inspect mesh quality."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "42e1c562",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%load_ext autoreload\n",
+    "%autoreload 2\n",
+    "\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import geopandas as gpd\n",
+    "from shapely.geometry import Point, Polygon, LineString, box\n",
+    "\n",
+    "from vorflow import ConceptualMesh, MeshGenerator, VoronoiTessellator\n",
+    "from vorflow.utils import calculate_mesh_quality, summarize_quality"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e2b415d8",
+   "metadata": {},
+   "source": [
+    "## 1. Define Geometry\n",
+    "We define a domain with a river, a fault, a well, and a refined zone."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "20f14eb4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Domain: 200x200 box\n",
+    "domain = box(0, 0, 200, 200)\n",
+    "# add a hole to the domain\n",
+    "hole = box(80, 80, 120, 120)\n",
+    "domain = Polygon(domain.exterior.coords, [hole.exterior.coords])\n",
+    "\n",
+    "\n",
+    "# Refined Zone: A smaller box inside\n",
+    "zone_poly = box(50, 50, 75, 75)\n",
+    "\n",
+    "# River: A sine wave\n",
+    "river_coords = [(x, 75 + 20 * np.sin(0.1 * x)) for x in range(-10, 220, 10)]\n",
+    "river_line = LineString(river_coords)\n",
+    "\n",
+    "# Fault: A straight line (Barrier)\n",
+    "fault_line = LineString([(100, 0), (100, 150)])\n",
+    "\n",
+    "# Well: A point\n",
+    "well_point1 = Point(25, 25)\n",
+    "well_point2 = Point(60,60)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "890b48eb",
+   "metadata": {},
+   "source": [
+    "## 2. Build Conceptual Mesh\n",
+    "Here we configure the mesh refinement.\n",
+    "- **resolution**: Target edge length.\n",
+    "- **dist_min**: Distance from the feature where resolution is constant.\n",
+    "- **dist_max**: Distance where resolution decays to background size."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "65365548",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cm = ConceptualMesh(crs=\"EPSG:3857\")\n",
+    "\n",
+    "# 1. Add Domain (Background)\n",
+    "cm.add_polygon(domain, zone_id=1, resolution=20.0)\n",
+    "\n",
+    "# 2. Add Refined Zone\n",
+    "# dist_max_out controls how fast the mesh grows outside this zone\n",
+    "cm.add_polygon(zone_poly,z_order=2,  zone_id=2, resolution=3.0, dist_max_out=20.0)\n",
+    "\n",
+    "\n",
+    "# 3. Add River (Standard Line)\n",
+    "cm.add_line(river_line, line_id=\"river\", resolution=2.0, dist_max=15.0)\n",
+    "\n",
+    "# 4. Add Fault (Barrier)\n",
+    "# is_barrier=True ensures the mesh edges align with this line\n",
+    "cm.add_line(fault_line, line_id=\"fault\", resolution=2.0, is_barrier=True, dist_max=15.0)\n",
+    "\n",
+    "# 5. Add Well\n",
+    "cm.add_point(well_point1, point_id=\"well-a\", resolution=1.0, dist_max=50.0)\n",
+    "cm.add_point(well_point2, point_id=\"well-b\", resolution=1.0, dist_max=20.0)\n",
+    "\n",
+    "# Process inputs\n",
+    "clean_polys, clean_lines, clean_pts = cm.generate()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "38c1334b",
+   "metadata": {},
+   "source": [
+    "## 3. Generate Triangular Mesh\n",
+    "We use Gmsh to generate the underlying triangulation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f1d756f6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mg = MeshGenerator(background_lc=20.0, verbosity=2)\n",
+    "mg.generate(clean_polys, clean_lines, clean_pts)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7860463e",
+   "metadata": {},
+   "source": [
+    "## 4. Voronoi Tessellation\n",
+    "Convert the triangular mesh to a Voronoi grid, clipped to the domain."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2a00aaba",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vt = VoronoiTessellator(mg, cm, clip_to_boundary=True)\n",
+    "grid = vt.generate()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "566b45a5",
+   "metadata": {},
+   "source": [
+    "## 5. Visualization\n",
+    "Plot the resulting grid along with the input features."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3db9769f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(12, 12))\n",
+    "ax.set_aspect('equal')\n",
+    "\n",
+    "# Plot Grid\n",
+    "grid.plot(ax=ax, edgecolor='black', linewidth=0.5, alpha=0.5, color='white')\n",
+    "\n",
+    "# Plot Features\n",
+    "gpd.GeoSeries([domain]).boundary.plot(ax=ax, color='black', linewidth=2, label='Domain')\n",
+    "gpd.GeoSeries([zone_poly]).boundary.plot(ax=ax, color='blue', linestyle='--', label='Refined Zone')\n",
+    "gpd.GeoSeries([river_line]).plot(ax=ax, color='cyan', linewidth=2, label='River')\n",
+    "gpd.GeoSeries([fault_line]).plot(ax=ax, color='red', linewidth=2, linestyle='-', label='Fault (Barrier)')\n",
+    "gpd.GeoSeries([well_point1]).plot(ax=ax, color='red', marker='*', markersize=100, label='Well')\n",
+    "gpd.GeoSeries([well_point2]).plot(ax=ax, color='red', marker='*', markersize=100,)\n",
+    "\n",
+    "plt.legend()\n",
+    "plt.title(\"Voronoi Grid with Refinement\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e8657d18",
+   "metadata": {},
+   "source": [
+    "## 6. Quality Analysis\n",
+    "Check the geometric quality of the grid cells."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0de14088",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "quality_gdf = calculate_mesh_quality(grid, calc_ortho=True)\n",
+    "summarize_quality(quality_gdf)\n",
+    "\n",
+    "# Plot Orthogonality Error\n",
+    "fig, ax = plt.subplots(figsize=(10, 8))\n",
+    "quality_gdf.plot(column='ortho_error', ax=ax, legend=True, cmap='Reds', vmin=0, vmax=10)\n",
+    "plt.title(\"Orthogonality Error (Degrees)\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8172e0e5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot Orthogonality Error\n",
+    "fig, ax = plt.subplots(figsize=(10, 8))\n",
+    "quality_gdf.plot(column='drift_ratio', ax=ax, legend=True, cmap='Reds', vmin=0, vmax=.4)\n",
+    "plt.title(\"Drift Ratio\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "857040ca",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "vorflow",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.14.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/src/vorflow/__pycache__/__init__.cpython-314.pyc b/src/vorflow/__pycache__/__init__.cpython-314.pyc
deleted file mode 100644
index 1427741..0000000
Binary files a/src/vorflow/__pycache__/__init__.cpython-314.pyc and /dev/null differ
diff --git a/src/vorflow/__pycache__/blueprint.cpython-314.pyc b/src/vorflow/__pycache__/blueprint.cpython-314.pyc
deleted file mode 100644
index fd71c65..0000000
Binary files a/src/vorflow/__pycache__/blueprint.cpython-314.pyc and /dev/null differ
diff --git a/src/vorflow/__pycache__/engine.cpython-314.pyc b/src/vorflow/__pycache__/engine.cpython-314.pyc
deleted file mode 100644
index f6c56a1..0000000
Binary files a/src/vorflow/__pycache__/engine.cpython-314.pyc and /dev/null differ
diff --git a/src/vorflow/__pycache__/tessellator.cpython-314.pyc b/src/vorflow/__pycache__/tessellator.cpython-314.pyc
deleted file mode 100644
index 34f48ac..0000000
Binary files a/src/vorflow/__pycache__/tessellator.cpython-314.pyc and /dev/null differ