Using Regions, can I model a reflecting wavefront?

You could use NDSolve to do the hard work:

region = Disk[];
sol = NDSolveValue[{D[u[t, x, y], {t, 2}] - 
     Laplacian[u[t, x, y], {x, y}] == 0, 
   DirichletCondition[u[t, x, y] == 0, True], 
   u[0, x, y] == 2*Exp[-125 ((x)^2 + (y - 0.5)^2)], 
   Derivative[1, 0, 0][u][0, x, y] == 0}, u, {t, 0, 2}, 
  Element[{x, y}, region]]

And then:

ListAnimate[
 Table[Plot3D[sol[t, x, y], Element[{x, y}, region], 
   PlotRange -> {-0.75, 2}, AspectRatio -> Automatic, Boxed -> False, 
   Axes -> None, PlotPoints -> 33], {t, 0, 2, 1/25}], 
 SaveDefinitions -> True]

To have an internal obstacle just change the region:

region = RegionDifference[Disk[], 
   Rectangle[{-1/3, -1/3}, {1/3, -1/4}]];
sol = NDSolveValue[{D[u[t, x, y], {t, 2}] - 
     Laplacian[u[t, x, y], {x, y}] == 0, 
   DirichletCondition[u[t, x, y] == 0, True], 
   u[0, x, y] == 2*Exp[-125 ((x)^2 + (y - 0.5)^2)], 
   Derivative[1, 0, 0][u][0, x, y] == 0}, u, {t, 0, 2}, 
  Element[{x, y}, region]]

Visualize:

ListAnimate[
 Table[
   Plot3D[sol[t, x, y], Element[{x, y}, region], 
    PlotRange -> {-0.75, 2}, AspectRatio -> Automatic, Boxed -> False,
     Axes -> None, PlotPoints -> 33], {t, 0, 2, 1/25}], 
 SaveDefinitions -> True]

Also, you can find much more information on the wave equation by looking at the Acoustics in the Time Domain tutorial in the documentation system under PDEModels/tutorial/AcousticsTimeDomain

I adapted @Kuba's approach from this answer to generate a quick and dirty particle tracer.

(* Create and Discretize Region *)
region = RegionDifference[Disk[], 
   Rectangle[{-1/3, -1/3}, {1/3, -1/4}]];
R2 = RegionBoundary@DiscretizeRegion@region;
rdf = RegionDistance[R2];
rnf = RegionNearest[R2];
(* Time Increment *)
dt = 0.001;
(* Collision Margin *)
margin = 1.05 dt;
r0 = 1000;
(* Starting Point for Emission *)
sp = {0, 0};
(* Conditional Particle Advancer *)
advance[r_, x_, v_, c_] := 
 Block[{xnew = x + dt v}, {rdf[xnew], xnew, v, c}] /; r > margin
advance[r_, x_, v_, c_] := 
 Block[{xnew = x , vnew = v, normal = Normalize[x - rnf[x]]},
   vnew = Normalize[v - 2 v.normal normal];
   xnew += dt vnew;
   {rdf[xnew], xnew, vnew, c + 1}] /; r <= margin

Now, we can run the simulation and create an animation at every 50 time steps.

nparticles = 1000;
ntimesteps = 5000;
tabres = Table[
   NestList[
    advance @@ # &, {rdf[sp], 
     sp, {Cos[2 Pi #], Sin[2 Pi #]} &@RandomReal[], 0}, 
    ntimesteps], {i, 1, nparticles}];
frames = Table[
   RegionPlot[R2, Epilog -> (Disk[#, 0.01] & /@ tabres[[All, i, 2]]), 
    AspectRatio -> Automatic], {i, 1, ntimesteps, 50}];
ListAnimate@frames

@user21's solution is very impressive. However, it isn't quite what I was looking for. This is because of the interaction between the waves. They are acting - well - like waves. This means we have a linear addition of the waves. This was what the original question forbid ;) . We want a single wavefront to come from the centre of the sphere and watch what happens as it moves around objects. Imagine it is a single photon - and doesn't act like a water wave.

Of course, if we're talking about single photons - a raytracing solution would work. I've implemented one (inspired from here), however, again - it isn't what the original question is asking for. We want a single wavefront which spreads....

But, anyway, this is my ray-tracing attempt

With 3 photons:

With 100 photons:

(* Line Intersection *) 


 LLI[vi_List] := 
 With[{x1 = vi[[1, 1]], y1 = vi[[1, 2]], x2 = vi[[2, 1]], 
   y2 = vi[[2, 2]], x3 = vi[[3, 1]], y3 = vi[[3, 2]], x4 = vi[[4, 1]],
    y4 = vi[[4, 
      2]]}, {-((-(x3 - x4) (x2 y1 - x1 y2) + (x1 - x2) (x4 y3 - 
           x3 y4))/((x3 - x4) (y1 - y2) + (-x1 + x2) (y3 - 
           y4))), (x4 (y1 - y2) y3 + x1 y2 y3 - x3 y1 y4 - x1 y2 y4 + 
      x3 y2 y4 + 
      x2 y1 (-y3 + y4))/(-(x3 - x4) (y1 - y2) + (x1 - x2) (y3 - y4))}]

(* Consider how we bounce *) 

bounce2[{p0_, d0_, i0_}] := 
 Module[{idxL, pL, validL, distL, i, p1, d1, bValid, dist, angleL, 
   angle}, idxL = 
   Position[Pi/2 < VectorAngle[d0, #] < Pi 3/2 Pi & /@ norm, True] // 
    Flatten;
  pL = Table[LLI[{p0, p0 + d0, ##}] & @@ edge[[j]], {j, idxL}];
  validL = 
   Table[! Or @@ (Greater[#, 
          1] & /@ (EuclideanDistance[#, pL[[i]]]/
            length[[idxL[[i]]]] & /@ edge[[idxL[[i]]]])), {i, 
     Length@idxL}];
  distL = EuclideanDistance[#, p0] & /@ pL;
  angleL = 
   Table[VectorAngle[norm[[idxL[[i]]]], pL[[i]] - p0], {i, 
     Length@idxL}];
  {i, p1, bValid, angle, dist} = 
   Select[Transpose@{idxL, pL, validL, angleL, 
        distL}, (#[[3]] && #[[4]] > Pi/2) &] // 
     MinimalBy[#, Last] & // #[[1]] &;
  d1 = (ReflectionTransform[RotationTransform[-Pi/2]@(-norm[[i]]), 
        p1]@p0 - p1) // Normalize;
  {p1, d1, i}]

(* Give our boundaries *) 

boundary1 = CirclePoints[2, 100];

edge1 = Table[
   RotateRight[boundary1, i][[;; 2]], {i, Length@boundary1}];
length1 = EuclideanDistance @@ # & /@ edge1;
norm1 = Normalize@(RotationTransform[Pi/2]@(#[[2]] - #[[1]])) & /@ 
   edge1;


boundary2 = {{-1, -0.2}, {1, -0.2}, {1, 0}, {-1, 0}};

edge2 = Table[
   RotateRight[boundary2, i][[;; 2]], {i, Length@boundary2}];
length2 = EuclideanDistance @@ # & /@ edge2;
norm2 = -Normalize@(RotationTransform[Pi/2]@(#[[2]] - #[[1]])) & /@ 
   edge2;

boundary = Join[boundary1, boundary2];
edge = Join[edge1, edge2];
length = Join[length1, length2];
norm = Join[norm1, norm2];


photons = 3; 
bounces = 100;
g = ConstantArray[{}, photons];

For[i = 1, i <= photons, i++, 
 p0 = {0, 0.1};
 d0 = {Cos@#, Sin@#} &@RandomReal[{0, 2 Pi}];
 r = NestList[bounce2, {p0, d0, 0}, bounces];
 p = r[[All, 1]];
 g[[i]] = 
  Table[Graphics[{FaceForm[LightBlue], EdgeForm[], Gray, 
     Line@p[[;; j]], Darker@Gray, Point@p[[;; j]], Red, 
     Point@p[[1]]}], {j, 2, Length@r}];
 ]

surface = 
 Graphics[{{FaceForm[LightBlue], Polygon@boundary1}, FaceForm[White], 
   Polygon@boundary2}]
animate = Table[Show[surface, g[[;; , {i}]]], {i, 1, bounces}];

ListAnimate[animate]

This isn't a complete solution as I really looking for the propagation of the circles around the sphere.