Crop a Voronoi diagram and get a proper MeshRegion

First Try (doesn't work due to a bug)

That's very frustrating. I thought that it might be a good idea to discretize the disk first and compute the intersections afterwards. That improves the timing by a factor of 20 but it gets some of the intersections wrong.

disk = BoundaryDiscretizeRegion[Disk[], MaxCellMeasure -> 0.01];
cells = MeshPrimitives[voronoi, 2];
pieces2 = BoundaryDiscretizeRegion@RegionIntersection[disk, #] & /@ cells; // AbsoluteTiming // First
Show[
 Cases[pieces2, _BoundaryMeshRegion],
 Graphics[{Red,
   Extract[cells, 
    Position[pieces2, Except[_BoundaryMeshRegion], 1, Heads -> False]]
   }]
 ]

1.45424

This is the result by Mmm 11.3. The intersections of the red cells with the circle could not be computed (probably a bug?):

Second Try

Not as simple anymore but quite efficient and maybe sufficiently general.

Let's say R is a MeshRegion and B is BoundaryMeshRegion. We want to crop R to be contained in B.

The strategy: Cells of R are divided into three groups: internalcells which are entirely contained in B, externalcells which do not intersect with B, and the bndcells which intersect the boundary. The internalcells can just be copied, the externalcells get discarded and the bndcells are first discretized by BoundaryDiscretizeRegion and then intersected with B. In the end, everythng is merged by the method that you fortunately mentioned. (I completely forgot about that post). Note that also IGMeshCellAdjacencyMatrix comes in handy here, in counting the number of vertices of each cell that lie within B. Note however that this is guaranteed to work correctly only if each cell of R is also convex (a mild condition) and if B is also convex (a more restrictive condition).

Needs["IGraphM`"];
cropRegion[R_MeshRegion, B_BoundaryMeshRegion] := 
 Module[{f, pts, cells, booleans, A, intersec, internalcells, 
   externalcells, bndcells, bndpieces, ϵ, plist, newpts, nf, 
   cellprimitives, celldata},
  f = RegionMember[B];
  pts = MeshCoordinates[R];
  cellprimitives = MeshPrimitives[R, 2];
  booleans = Map[Boole@*f, pts];
  
  A = IGMeshCellAdjacencyMatrix[R, 2, 0];
  intersec = (A.booleans);
  
  internalcells = Flatten[Position[Subtract[intersec, Total[A, {2}]], 0]];
  externalcells = Flatten[Position[intersec, 0]];
  bndcells = Complement[Range[Length[A]], internalcells, externalcells];
  
  bndpieces = ParallelMap[
     RegionIntersection[B, With[{p = #[[1]]}, BoundaryMesh@MeshRegion[p, Head[#][Range[Length[p]]]]]] &, 
     cellprimitives[[bndcells]]
     ];
  
  celldata = Join[
    Join @@ (MeshPrimitives[#, 2] & /@ bndpieces)[[All, All, 1]],
    cellprimitives[[internalcells, 1]]
    ];
  ϵ = 0.5 Min[
     PropertyValue[{R, 1}, MeshCellMeasure],
     PropertyValue[{#, 1}, MeshCellMeasure] & /@ bndpieces
     ]; 
  pts = Flatten[celldata, 1]; 
  plist = Sort[DeleteDuplicates[
     Compile[{{idx, _Integer, 1}}, 
       Min[idx], 
       RuntimeAttributes -> {Listable}, 
       Parallelization -> True
       ][Nearest[pts -> Automatic, pts, {∞, ϵ}]]
     ]]; 
  newpts = pts[[plist]]; nf = Nearest[newpts -> Automatic]; 
  cells = Internal`PartitionRagged[Flatten[nf[pts]], Length /@ celldata];
  MeshRegion[newpts, Polygon[cells]]
  ];

Usage example and timing

voronoi = VoronoiMesh@RandomPoint[Disk[], 10000];
r = 1;
disk = BoundaryDiscretizeRegion[Disk[{0, 0}, r], MaxCellMeasure -> 0.01];
S = cropRegion[voronoi, disk]; // RepeatedTiming // First
S

0.416

The original Voronoi mesh with 100 cells gets processed in about 0.04 seconds. That's about 700 times as fast as what we have started from.

There is still some optimization potential since we do not really need to search the interior vertices of the interior region for duplicates...

Warning

If B is not convex, then some less pleasant things can happen. For example, using this star fish

B = BoundaryMeshRegion[
  Map[t \[Function] (2 + Cos[5 t])/3 {Cos[t], Sin[t]}, Most@Subdivide[0., 2. Pi, 2000]],
  Line[Partition[Range[2000], 2, 1, 1]]
  ]

can lead to

SeedRandom[17];
cropRegion[VoronoiMesh@RandomPoint[Disk[], 100], B]

But this becomes more and more negligible, the smaller the cells in R are:

cropRegion[VoronoiMesh@RandomPoint[Disk[], 10000], B]

pieces2 = BoundaryDiscretizeRegion /@ 
  (Graphics`PolygonUtils`PolygonIntersection[#, Polygon[CirclePoints[100]]] &/@ 
    MeshPrimitives[voronoi, 2]); // AbsoluteTiming // First

1.09499

(vs 33.0193 for pieces = ...)

Show[pieces2]

With 1000 points it takes 11.3679 seconds to get pieces2, and Show[pieces2] gives:

Here's another implementation based on the idea of grabbing cells fully contained in the disk and intersecting the ones with partial overlap.

maskedMeshRegion[mr_?MeshRegionQ, mask_?ConstantRegionQ, δ_:0.01] :=
  Block[{prims, inside, frontier, bmr, fpieces, pts, cells},
    prims = MeshPrimitives[mr, 2];

    inside = Select[prims, RegionWithin[mask, #]&];

    frontier = Select[Complement[prims, inside], !RegionDisjoint[mask, #]&];
    bmr = BoundaryDiscretizeRegion[mask, MaxCellMeasure -> δ];
    fpieces = RegionIntersection[bmr, 
      BoundaryMeshRegion[#, Line[Append[Range[Length[#]], 1]]]]& /@ frontier[[All, 1]];

    pts = Union[Join @@ MeshCoordinates /@ fpieces, Join @@ inside[[All, 1]]];
    cells = Join[
      Join @@ (MeshPrimitives[#, 2]& /@ fpieces), 
      inside
    ] /. Dispatch[Thread[pts -> Range[Length[pts]]]];

    MeshRegion[pts, cells]
]

A test:

voronoi = VoronoiMesh@RandomPoint[Disk[], 10000];
maskedMeshRegion[voronoi, Disk[]] // AbsoluteTiming