diff --git a/kwant/builder.py b/kwant/builder.py
index 2e1514b8dde403973fe519c6657e3c4608c15a04..f5c9f023a88b90427d9bacc6c7e2e0fc982bcf4e 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -1263,8 +1263,17 @@ class Builder:
     def fill(self, template, shape, start, *, max_sites=10**7):
         """Populate builder using another one as a template.
 
-        Sites and hoppings that already exist in the target builder are never
-        overwritten.
+        Starting from one or multiple sites, traverse the graph of the template
+        builder and copy sites and hoppings to the target builder.  The
+        traversal stops at sites that are already present in the target and on
+        sites that are not inside the provided shape.
+
+        This function takes into account translational symmetry.  As such,
+        typically the template will have a higher symmetry than the target.
+
+        Newly added sites are connected by hoppings to sites that were already
+        present.  This facilitates construction of a system by a series of
+        calls to 'fill'.
 
         Parameters
         ----------
@@ -1380,21 +1389,21 @@ class Builder:
                         if (head_fd not in old_active
                             and head_fd not in new_active):
                             # The 'head' site has not been filled yet.
-                            if not shape(head_fd):
-                                continue
-
-                            if head_fd not in H:
-                                # Fill 'head' site.
-                                new_active.add(head_fd)
-                                H.setdefault(head_fd, [head_fd, None])
-                            else:
-                                # The 'head' site exists and won't be visited:
-                                # fill the incoming edge as well to balance the
-                                # hopping.
+                            if head_fd in H:
+                                # The 'head' site exists.  (It doesn't matter
+                                # whether it's in the shape or not.)  Fill the
+                                # incoming edge as well to balance the hopping.
                                 other_value = template._get_edge(
                                     *templ_sym.to_fd(head, tail))
                                 self._set_edge(*to_fd(head, tail)
                                                + (other_value,))
+                            else:
+                                if not shape(head_fd):
+                                    # There is no site at 'head' and it's
+                                    # outside the shape.
+                                    continue
+                                new_active.add(head_fd)
+                                H.setdefault(head_fd, [head_fd, None])
 
                         # Fill the outgoing edge.
                         if head in old_heads:
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index 8eb2e58b9ced9e9f32038e3776dacbe737eef4ef..b31c0dc09cf18cf96099aac1e69e33ebaf23b0c3 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -790,6 +790,42 @@ def test_fill():
     assert to_target_fd(new_sites[0]) == to_target_fd(lat(-1))
 
 
+def test_fill_sticky():
+    """Test that adjacent regions are properly interconnected when filled
+    separately.
+    """
+    # Generate model template.
+    lat = kwant.lattice.kagome()
+    template = kwant.Builder(kwant.TranslationalSymmetry(
+        lat.vec((1, 0)), lat.vec((0, 1))))
+    for i, sl in enumerate(lat.sublattices):
+        template[sl(0, 0)] = i
+    for i in range(1, 3):
+        for j, hop in enumerate(template.expand(lat.neighbors(i))):
+            template[hop] = j * 1j
+
+    def disk(site):
+        pos = site.pos
+        return ta.dot(pos, pos) < 13
+
+    def halfplane(site):
+        return ta.dot(site.pos - (-1, 1), (-0.9, 0.63)) > 0
+
+    # Fill in one go.
+    syst0 = kwant.Builder()
+    syst0.fill(template, disk, (0, 0))
+
+    # Fill in two stages.
+    syst1 = kwant.Builder()
+    syst1.fill(template, lambda s: disk(s) and halfplane(s), (-2, 1))
+    syst1.fill(template, lambda s: disk(s) and not halfplane(s), (0, 0))
+
+    # Verify that both results are identical.
+    assert set(syst0.site_value_pairs()) == set(syst1.site_value_pairs())
+    assert (set(syst0.hopping_value_pairs())
+            == set(syst1.hopping_value_pairs()))
+
+
 def test_attach_lead():
     fam = builder.SimpleSiteFamily()
     fam_noncommensurate = builder.SimpleSiteFamily(name='other')