diff --git a/src/vorflow/blueprint.py b/src/vorflow/blueprint.py
index bb869a9..bb524b6 100644
--- a/src/vorflow/blueprint.py
+++ b/src/vorflow/blueprint.py
@@ -4,6 +4,11 @@
 from shapely.geometry import Polygon, LineString, Point, box, MultiPolygon
 from shapely.ops import unary_union, snap, linemerge
 from shapely.validation import make_valid
+from shapely.strtree import STRtree
+
+# Constants for geometry simplification and reporting
+SIGNIFICANT_REDUCTION_PCT = 1.0
+DEFAULT_TOLERANCE = 1e-3
 
 class ConceptualMesh:
     def __init__(self, crs="EPSG:4326"):
@@ -29,7 +34,20 @@ def __init__(self, crs="EPSG:4326"):
         self.clean_lines = gpd.GeoDataFrame()
         self.clean_points = gpd.GeoDataFrame()
 
-    def add_polygon(self, geometry, zone_id, resolution=None, z_order=0, mesh_refinement=True, dist_min=None, dist_max=None, dist_max_in=None, dist_max_out=None, border_density=None):
+    def add_polygon(
+        self,
+        geometry,
+        zone_id,
+        resolution=None,
+        z_order=0,
+        mesh_refinement=True,
+        dist_min=None,
+        dist_max=None,
+        dist_max_in=None,
+        dist_max_out=None,
+        densify=None,
+        simplify_tolerance=None,
+    ):
         """
         Adds a polygon feature, such as a model boundary or a refinement zone.
 
@@ -49,29 +67,52 @@ def add_polygon(self, geometry, zone_id, resolution=None, z_order=0, mesh_refine
                 transitions from the boundary resolution to the internal resolution.
             dist_max_out (float, optional): Distance outside the polygon over which the
                 mesh transitions to the background resolution.
-            border_density (float, optional): If set, densifies the polygon's boundary
-                by adding vertices, ensuring no segment is longer than this value.
+            densify (float|bool|None, optional): Controls polygon boundary densification:
+                - If False, disables densification.
+                - If True, densifies using `resolution` (lc). Requires `resolution` to be set.
+                - If a float, densifies so no boundary segment is longer than this value.
+                Raises ValueError if non-positive when specified as a float.
+            simplify_tolerance (float|int|None, optional): If a number > 0, applies Douglas-Peucker
+                simplification with this tolerance. If None or 0, no simplification is applied.
+                Raises ValueError if negative. Boolean values are not supported.
         """
         if not geometry.is_valid:
             geometry = make_valid(geometry)
-            
+
+        if isinstance(simplify_tolerance, bool):
+            raise ValueError(
+                "simplify_tolerance must be a non-negative number (or None/0 to disable). "
+                "Boolean values are not supported."
+            )
+        if isinstance(simplify_tolerance, (int, float)) and simplify_tolerance < 0:
+            raise ValueError(f"simplify_tolerance must be non-negative. Got {simplify_tolerance}.")
+
+        if isinstance(densify, (int, float)) and not isinstance(densify, bool) and densify <= 0:
+            raise ValueError(f"densify must be positive when specified as a float. Got {densify}.")
+        if densify is True and (resolution is None or resolution <= 0):
+            raise ValueError("densify=True for polygons requires a positive `resolution` (lc).")
+
         # For backward compatibility, allow 'dist_max' to function as 'dist_max_out'.
         if dist_max is not None and dist_max_out is None:
             dist_max_out = dist_max
 
-        self.raw_polygons.append({
-            'geometry': geometry,
-            'zone_id': zone_id,
-            'lc': resolution,
-            'z_order': z_order,
-            'refine': mesh_refinement,
-            'dist_min': dist_min,
-            'dist_max_in': dist_max_in,
-            'dist_max_out': dist_max_out,
-            'border_density': border_density
-        })
+        self.raw_polygons.append(
+            {
+                "geometry": geometry,
+                "zone_id": zone_id,
+                "lc": resolution,
+                "z_order": z_order,
+                "refine": mesh_refinement,
+                "dist_min": dist_min,
+                "dist_max_in": dist_max_in,
+                "dist_max_out": dist_max_out,
+                "densify": densify,
+                "simplify_tolerance": simplify_tolerance,
+            }
+        )
 
-    def add_line(self, geometry, line_id, resolution, snap_to_polygons=True, is_barrier=False, dist_min=None, dist_max=None, straddle_width=None):
+    def add_line(self, geometry, line_id, resolution, snap_to_polygons=True, is_barrier=False,
+                  dist_min=None, dist_max=None, straddle_width=None, densify=True, simplify_tolerance=None):
         """
         Adds a line feature, such as a river, fault, or other linear boundary.
 
@@ -89,10 +130,29 @@ def add_line(self, geometry, line_id, resolution, snap_to_polygons=True, is_barr
                 transitions to the background resolution.
             straddle_width (float, optional): If set, forces Voronoi cell edges to align
                 perfectly with the line by creating a "virtual straddle" of mesh nodes.
+            densify (float or bool, optional): Controls line densification:
+                - If False, disables densification.
+                - If True, densifies the line using the `resolution` value.
+                - If a float, densifies the line so that no segment is longer than this value.
+                Raises ValueError if negative or zero.
+            simplify_tolerance (float|int|None, optional): If a number > 0, simplifies the line with
+                this tolerance using Douglas-Peucker algorithm. If None or 0, no simplification is applied.
+                Raises ValueError if negative. Boolean values are not supported.
         """
         if not geometry.is_valid:
             geometry = make_valid(geometry)
-            
+
+        if isinstance(simplify_tolerance, bool):
+            raise ValueError(
+                "simplify_tolerance must be a non-negative number (or None/0 to disable). "
+                "Boolean values are not supported."
+            )
+        if isinstance(simplify_tolerance, (int, float)) and simplify_tolerance < 0:
+            raise ValueError(f"simplify_tolerance must be non-negative. Got {simplify_tolerance}.")
+        
+        if isinstance(densify, (int, float)) and not isinstance(densify, bool) and densify <= 0:
+            raise ValueError(f"densify must be positive when specified as a float. Got {densify}.")
+        
         self.raw_lines.append({
             'geometry': geometry,
             'line_id': line_id,
@@ -100,10 +160,12 @@ def add_line(self, geometry, line_id, resolution, snap_to_polygons=True, is_barr
             'is_barrier': is_barrier,
             'dist_min': dist_min,
             'dist_max': dist_max,
-            'straddle_width': straddle_width
+            'straddle_width': straddle_width,
+            'densify': densify,
+            'simplify_tolerance': simplify_tolerance
         })
 
-    def add_point(self, geometry, point_id, resolution, dist_min=None, dist_max=None):
+    def add_point(self, geometry, point_id, resolution, dist_min=None, dist_max=None, simplify_tolerance=None):
         """
         Adds a point feature, such as a well or an observation point.
 
@@ -115,14 +177,130 @@ def add_point(self, geometry, point_id, resolution, dist_min=None, dist_max=None
                 held constant at the point's resolution.
             dist_max (float, optional): Distance from the point over which the mesh
                 transitions to the background resolution.
+            simplify_tolerance (float|int|None, optional): If a number > 0, merges points that are closer
+                than this tolerance. If None or 0, no merging is applied. Raises ValueError if negative.
+                Boolean values are not supported.
         """
+        if isinstance(simplify_tolerance, bool):
+            raise ValueError(
+                "simplify_tolerance must be a non-negative number (or None/0 to disable). "
+                "Boolean values are not supported."
+            )
+        if isinstance(simplify_tolerance, (int, float)) and simplify_tolerance < 0:
+            raise ValueError(f"simplify_tolerance must be non-negative. Got {simplify_tolerance}.")
         self.raw_points.append({
             'geometry': geometry,
             'point_id': point_id,
             'lc': resolution,
             'dist_min': dist_min,
-            'dist_max': dist_max
+            'dist_max': dist_max,
+            'simplify_tolerance': simplify_tolerance
         })
+    def _apply_simplification(self):
+        """
+        Applies geometry simplification to raw polygons, lines, and points
+        based on their specified tolerances to reduce geometric complexity.
+
+        For polygons and lines the Douglas-Peucker algorithm is used.
+        For points, this method performs deduplication: points that are within
+        a specified tolerance of each other are merged, and only the point with the
+        finest (smallest) resolution is kept.
+        """
+        # Simplify Polygons
+        for i, poly_data in enumerate(self.raw_polygons):
+            tol = poly_data.get('simplify_tolerance')
+            if isinstance(tol, bool):
+                raise ValueError(
+                    "simplify_tolerance must be a non-negative number (or None/0 to disable). "
+                    "Boolean values are not supported."
+                )
+
+            if tol is not None and tol > 0:
+                org_area = poly_data['geometry'].area
+                simplified_geom = poly_data['geometry'].simplify(tol, preserve_topology=True)
+                self.raw_polygons[i]['geometry'] = simplified_geom
+                new_area = simplified_geom.area
+                if new_area < org_area and org_area > 0:
+                    reduction_pct = 100 * (org_area - new_area) / org_area
+                    if reduction_pct > SIGNIFICANT_REDUCTION_PCT:
+                        print(
+                            f"Simplified polygon (zone_id={poly_data['zone_id']}) "
+                            f"reduced area by {reduction_pct:.2f}% using tolerance {tol}."
+                        )
+
+        # Simplify Lines
+        for i, line_data in enumerate(self.raw_lines):
+            tol = line_data.get('simplify_tolerance')
+            if isinstance(tol, bool):
+                raise ValueError(
+                    "simplify_tolerance must be a non-negative number (or None/0 to disable). "
+                    "Boolean values are not supported."
+                )
+
+            if tol is not None and tol > 0:
+                org_length = line_data['geometry'].length
+                simplified_geom = line_data['geometry'].simplify(tol, preserve_topology=True)
+                self.raw_lines[i]['geometry'] = simplified_geom
+                new_length = simplified_geom.length
+                if new_length < org_length and org_length > 0:
+                    reduction_pct = 100 * (org_length - new_length) / org_length
+                    if reduction_pct > SIGNIFICANT_REDUCTION_PCT:
+                        print(
+                            f"Simplified line (line_id={line_data['line_id']}) "
+                            f"reduced length by {reduction_pct:.2f}% using tolerance {tol}."
+                        )
+
+        # Merge points that are very close to each other (deduplication)
+        if self.raw_points:
+            # Let's sort by resolution first
+            sorted_points = sorted(
+                self.raw_points,
+                key=lambda x: x['lc'] if x['lc'] is not None else float('inf')
+            )
+            final_points = []
+            geoms = [p['geometry'] for p in sorted_points]
+            tree = STRtree(geoms)
+            kept_indices = set()
+
+            for i, point_data in enumerate(sorted_points):
+                current_geom = point_data['geometry']
+                tol = point_data.get('simplify_tolerance')
+
+                if isinstance(tol, bool):
+                    raise ValueError(
+                        "simplify_tolerance must be a non-negative number (or None/0 to disable). "
+                        "Boolean values are not supported."
+                    )
+
+                # None or <=0 => no merging for this point (keep as-is)
+                if tol is None or tol <= 0:
+                    final_points.append(point_data)
+                    kept_indices.add(i)
+                    continue
+
+                is_merged = False
+
+                # Query tree for potential neighbors
+                # tree.query returns indices of geometries that intersect the buffer
+                search_area = current_geom.buffer(tol)
+                candidate_indices = tree.query(search_area)
+
+                for candidate_idx in candidate_indices:
+                    if candidate_idx in kept_indices:
+                        if geoms[candidate_idx].distance(current_geom) < tol:
+                            is_merged = True
+                            break
+
+                if not is_merged:
+                    final_points.append(point_data)
+                    kept_indices.add(i)
+
+            if len(self.raw_points) != len(final_points):
+                print(
+                    f"Simplification merged {len(self.raw_points) - len(final_points)} "
+                    f"points out of {len(self.raw_points)}"
+                )
+            self.raw_points = final_points
 
     def _resolve_overlaps(self):
         """
@@ -131,6 +309,24 @@ def _resolve_overlaps(self):
         lower ones.
         """
         # Sort polygons by priority, with the highest z_order processed first.
+        if not self.raw_polygons:
+            # If this is empty, just create an empty GeoDataFrame.
+            self.clean_polygons = gpd.GeoDataFrame(
+                columns=[
+                    "geometry",
+                    "zone_id",
+                    "lc",
+                    "z_order",
+                    "refine",
+                    "dist_min",
+                    "dist_max_in",
+                    "dist_max_out",
+                    "densify",
+                    "simplify_tolerance",
+                ],
+                crs=self.crs,
+            )
+            return
         df = pd.DataFrame(self.raw_polygons)
         df = df.sort_values(by='z_order', ascending=False)
         
@@ -180,7 +376,7 @@ def _resolve_overlaps(self):
 
         self.clean_polygons = gpd.GeoDataFrame(final_features, crs=self.crs)
 
-    def _enforce_connectivity(self, tolerance=1e-3):
+    def _enforce_connectivity(self, tolerance=DEFAULT_TOLERANCE):
         """
         Snaps features together to ensure they are topologically connected before
         being passed to the mesher. This is crucial for Gmsh to correctly
@@ -231,6 +427,9 @@ def generate(self):
         ensures topological connectivity, and prepares clean GeoDataFrames
         for the mesher.
         """
+        print("Applying optional geometry simplification...")
+        self._apply_simplification()
+
         print("Resolving polygon overlaps...")
         self._resolve_overlaps()
         
@@ -310,15 +509,43 @@ def densify_line(line, max_segment_length):
 
     def _apply_densification(self):
         """Applies densification to the clean polygon and line features."""
-        # Densify polygon boundaries where a 'border_density' is specified.
+        # Densify polygon boundaries based on `densify`.
         if not self.clean_polygons.empty:
-            self.clean_polygons['geometry'] = self.clean_polygons.apply(
-                lambda row: self._densify_geometry(row['geometry'], row['border_density']) 
-                if pd.notna(row.get('border_density')) else row['geometry'], axis=1
-            )
+
+            def get_poly_resolution(row):
+                d = row.get("densify")
+                if d is False or pd.isna(d):
+                    return None
+                if d is True:
+                    return row.get("lc")
+                if isinstance(d, (int, float)) and not isinstance(d, bool) and d > 0:
+                    return d
+                return None
+
+            def _poly_densify(row):
+                res = get_poly_resolution(row)
+                return self._densify_geometry(row["geometry"], res) if res is not None else row["geometry"]
+
+            self.clean_polygons["geometry"] = self.clean_polygons.apply(_poly_densify, axis=1)
 
         # Densify lines based on their target resolution ('lc').
         if not self.clean_lines.empty:
-            self.clean_lines['geometry'] = self.clean_lines.apply(
-                lambda row: self._densify_geometry(row['geometry'], row['lc']), axis=1
-            )
\ No newline at end of file
+            # Helper to determine the target resolution for a line row
+            def get_line_resolution(row):
+                d = row.get('densify')
+                
+                # 1. Explicitly disabled (densify=False)
+                if d is False:
+                    return None
+                
+                # 2. Explicit custom resolution (e.g., densify=5.0)
+                if isinstance(d, (int, float)) and not isinstance(d, bool) and d > 0:
+                    return d
+                
+                # 3. Default behavior (True or None): use the mesh resolution (lc)
+                return row.get('lc')
+            def _line_densify(row):
+                res = get_line_resolution(row)
+                return self._densify_geometry(row['geometry'], res) if res is not None else row['geometry']
+
+            self.clean_lines['geometry'] = self.clean_lines.apply(_line_densify, axis=1)
\ No newline at end of file
diff --git a/src/vorflow/engine.py b/src/vorflow/engine.py
index 942cde2..ed66d4d 100644
--- a/src/vorflow/engine.py
+++ b/src/vorflow/engine.py
@@ -494,10 +494,14 @@ def add_refinement(entity_dim, entity_tags, size_target, dist_min, dist_max, siz
         for idx, row in polygons_gdf.iterrows():
             if idx in gmsh_map['surfaces']:
                 tags = extract_tags(gmsh_map['surfaces'][idx])
-                
+
                 target_lc = get_row_param(row, 'lc', global_max_lc)
-                border_dens = get_row_param(row, 'border_density', target_lc)
-                boundary_lc = min(target_lc, border_dens)
+
+                densify_val = row.get("densify", None)
+                if isinstance(densify_val, (int, float)) and not isinstance(densify_val, bool) and densify_val > 0:
+                    boundary_lc = min(target_lc, float(densify_val))
+                else:
+                    boundary_lc = target_lc
 
                 if self.verbosity > 1:
                     print(f"Poly {idx}: Target={target_lc}, Border={boundary_lc}, Global={global_max_lc}")
diff --git a/tests/test_conceptual_mesh.py b/tests/test_conceptual_mesh.py
index 6c8388b..2ced0bb 100644
--- a/tests/test_conceptual_mesh.py
+++ b/tests/test_conceptual_mesh.py
@@ -49,3 +49,112 @@ def test_lines_and_points_snap_to_polygons():
     tolerance = 1e-3
     assert snapped_line.distance(boundary) <= tolerance
     assert snapped_point.distance(boundary) <= tolerance
+
+def test_polygon_simplification():
+    """Test that polygons are simplified when tolerance is provided."""
+    cm = ConceptualMesh()
+    
+    # Create a "noisy" square with a tiny bump on the top edge
+    # (0,1) -> (0.5, 1.001) -> (1,1)
+    poly = Polygon([(0, 0), (1, 0), (1, 1), (0.5, 1.001), (0, 1)])
+    
+    # Add with a tolerance larger than the noise (0.001)
+    cm.add_polygon(poly, zone_id=1, simplify_tolerance=0.01)
+    
+    clean_polys, _, _ = cm.generate()
+    
+    simplified_geom = clean_polys.iloc[0].geometry
+    
+    # The original polygon has 5 vertices + closing = 6 points in exterior ring
+    # The simplified one should remove the bump, leaving 4 corners + closing = 5 points
+    assert len(simplified_geom.exterior.coords) == 5
+    assert len(simplified_geom.exterior.coords) < len(poly.exterior.coords)
+
+def test_line_simplification():
+    """Test that lines are simplified when tolerance is provided."""
+    cm = ConceptualMesh()
+    # Noisy line: straight but with a midpoint slightly off
+    line = LineString([(0, 0), (0.5, 0.001), (1, 0)])
+    
+    cm.add_line(line, line_id="noisy_line", resolution=0.1, simplify_tolerance=0.01, densify=False)
+    
+    _, clean_lines, _ = cm.generate()
+    
+    simplified_line = clean_lines.iloc[0].geometry
+    # Should be simplified to just start and end points
+    assert len(simplified_line.coords) == 2
+
+def test_point_deduplication():
+    """Test that close points are merged and the finest resolution is kept."""
+    cm = ConceptualMesh()
+    p1 = Point(0, 0)
+    p2 = Point(0.0001, 0) # Very close to p1
+    
+    # Case 1: No simplification (default) -> Should keep both
+    cm.add_point(p1, "p1", resolution=1.0)
+    cm.add_point(p2, "p2", resolution=0.5) 
+    
+    _, _, clean_points = cm.generate()
+    assert len(clean_points) == 2
+    
+    # Case 2: With simplification -> Should merge
+    cm2 = ConceptualMesh()
+    # p2 has finer resolution (0.5), so it should be the one kept
+    cm2.add_point(p1, "p1", resolution=1.0, simplify_tolerance=0.01)
+    cm2.add_point(p2, "p2", resolution=0.5, simplify_tolerance=0.01)
+    
+    _, _, clean_points_merged = cm2.generate()
+    
+    assert len(clean_points_merged) == 1
+    
+    # Verify we kept the point with the finer resolution (0.5)
+    kept_point = clean_points_merged.iloc[0]
+    assert kept_point['lc'] == 0.5
+    assert kept_point['point_id'] == "p2"
+
+def test_line_densification_options():
+    """Test the three modes of line densification: False, True, and float."""
+    cm = ConceptualMesh()
+    # A line of length 10
+    line = LineString([(0, 0), (10, 0)])
+    
+    # 1. densify=False: Should NOT add vertices
+    cm.add_line(line, "no_densify", resolution=1.0, densify=False)
+    
+    # 2. densify=True (default): Should use resolution (1.0) -> ~10 segments
+    cm.add_line(line, "default_densify", resolution=1.0, densify=True)
+    
+    # 3. densify=5.0: Should use custom spacing (5.0) -> ~2 segments
+    cm.add_line(line, "custom_densify", resolution=1.0, densify=5.0)
+    
+    _, clean_lines, _ = cm.generate()
+    
+    # Check 1: No densification
+    l1 = clean_lines[clean_lines['line_id'] == "no_densify"].iloc[0].geometry
+    assert len(l1.coords) == 2 # Just start and end
+    
+    # Check 2: Default densification (lc=1.0)
+    l2 = clean_lines[clean_lines['line_id'] == "default_densify"].iloc[0].geometry
+    # Should have 11 points (10 segments)
+    assert len(l2.coords) == 11 
+    
+    # Check 3: Custom densification (val=5.0)
+    l3 = clean_lines[clean_lines['line_id'] == "custom_densify"].iloc[0].geometry
+    # Should have roughly 3 points (2 segments)
+    assert len(l3.coords) == 3
+
+@pytest.mark.parametrize("bool_tol", [True, False])
+def test_simplify_tolerance_bool_is_rejected(bool_tol):
+    cm = ConceptualMesh()
+
+    poly = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
+    with pytest.raises(ValueError):
+        cm.add_polygon(poly, zone_id=1, simplify_tolerance=bool_tol)
+
+    line = LineString([(0, 0), (1, 0)])
+    with pytest.raises(ValueError):
+        cm.add_line(line, line_id="l1", resolution=0.1, simplify_tolerance=bool_tol, densify=False)
+
+    pt = Point(0, 0)
+    with pytest.raises(ValueError):
+        cm.add_point(pt, point_id="p1", resolution=0.1, simplify_tolerance=bool_tol)
diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py
index a814166..e477a26 100644
--- a/tests/test_pipeline.py
+++ b/tests/test_pipeline.py
@@ -55,7 +55,7 @@ def __init__(self):
 def _build_simple_conceptual_mesh():
     cm = ConceptualMesh(crs="EPSG:3857")
     square = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
-    cm.add_polygon(square, zone_id=99, border_density=0.5)
+    cm.add_polygon(square, zone_id=99, densify=0.5)
     return cm