ContourPlot of a discontinuous function

Question

rainbow = {"Rainbow", {1.7, 3.4}};
color = ColorData[rainbow];
cf = Function[{z}, Directive[Opacity[0.8], ColorData[rainbow][z]]];
plot = {AspectRatio -> 1,
ImageSize -> {400, 400},
PlotRange -> {{-1, 1}, {-1, 1}},
FrameLabel -> {x, y},
LabelStyle -> Black,
FrameStyle -> Directive[Thickness[0.002], FontSize -> 14],
Mesh -> None, MeshStyle -> Gray,
InterpolationOrder -> 0,
ColorFunctionScaling -> False, ColorFunction -> cf};

g[x_, y_, z_] := Piecewise[{{(x - 1/2)^2 + y^2 + z^2 - 1/16, (x - 1/2)^2 + y^2 - 1/16 <= 0 && z >= 0}}, z];
(*t is a function of x1 and x2, and it is a discontinuous function as shown in the Figure 1*)
t[x1_, x2_] := Module[
            {center, line, surface, intersection, interPoint, d, data},
            center = {0, 0, 1};
            line = Line[{{x1, x2, 0}, center}];
            surface = ImplicitRegion[g[x, y, z] == 0., {x, y, z}];
            intersection = MeshCoordinates@DiscretizeRegion@RegionIntersection[line, surface];
            interPoint = Nearest[intersection, center][[1]];
            d = EuclideanDistance[center, interPoint];
            time = 2 d;
            time
            ];
data = Table[t[x, y], {y, -1, 1.0, 0.05}, {x, -1, 1.0, 0.05}];
ListDensityPlot[data, DataRange -> {{-1, 1}, {-1, 1}}, Evaluate@plot]

You don't need to read Module carefully, t[x,y] is only a function of x and y and it is discontinuous at the region as shown in Figure.

Then I want to plot some contourplot lines. For example:

ContourPlot[t[x, y] == 2.3, {x, -1, 1}, {y, -1, 1}, PlotRange -> {{-1, 1}, {-1, 1}}, PlotPoints -> 20, MaxRecursion -> 1]

and it looks like Figure 2(a). But what I want is Figure2(b), because of the discontinuous region.

1. How to modify the code to plot Figure2(b)?
2. There is no problem when I run the ListDensityPlot , however, there is a warning when I run ContourPlot:

Nearest::near1: "MeshCoordinates[DiscretizeRegion[RegionIntersection[ImplicitRegion[Piecewise[{{<<2>>}},z]==0.,{x,y,z}],Line[{{x,y,0},{0,0,1}}]]]] is neither a list of real points nor a valid list of rules"

Why?

3. In the ListDensityPlot If the code is

   data = Table[t[x, y], {x, -1, 1.0, 0.05}, {y, -1, 1.0, 0.05}];

The plot is rotated 90 degress with what it should be. I have to change the code to

data = Table[t[x, y], {y, -1, 1.0, 0.05}, {x, -1, 1.0, 0.05}];

Why?

Problem 2 and 3 are not important. They didn't have an effect on the results, I am just curious about it.

This is my solution:

Show[ContourPlot[t[x, y] == 1.8, {x, 0.3, 0.7}, {y, -0.2, 0.2}, ContourStyle -> color[1.8], PlotPoints -> 20, MaxRecursion -> 0,  PlotRange -> {{-1, 1}, {-1, 1}}],
ContourPlot[t[x, y] == 1.9, {x, 0.2, 0.8}, {y, -0.3, 0.3},  ContourStyle -> color[1.9], PlotPoints -> 20, MaxRecursion -> 0],
ContourPlot[t[x, y] == 2.0, {x, 0.2, 0.8}, {y, -0.3, 0.3}, ContourStyle -> color[2.0], PlotPoints -> 20, MaxRecursion -> 1, Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.1, {x, -0.33, 0.85}, {y, -0.33, 0.33},  ContourStyle -> color[2.1], PlotPoints -> 20, MaxRecursion -> 1,  Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.2, {x, -0.46, 0.9}, {y, -0.46, 0.46},  ContourStyle -> color[2.2], PlotPoints -> 20, MaxRecursion -> 0,  Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.3, {x, -0.57, 0.9}, {y, -0.57, 0.57}, ContourStyle -> color[2.3], PlotPoints -> 20, MaxRecursion -> 0,  Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.4, {x, -0.67, 0.9}, {y, -0.67, 0.67},  ContourStyle -> color[2.4], PlotPoints -> 20, MaxRecursion -> 0, Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.5, {x, -0.76, 0.9}, {y, -0.76, 0.76}, ContourStyle -> color[2.5], PlotPoints -> 20, MaxRecursion -> 0,  Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.6, {x, -0.84, 0.84}, {y, -0.84, 0.84}, ContourStyle -> color[2.6], PlotPoints -> 20, MaxRecursion -> 0, Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}]]

I define the color=ColorData[{"Rainbow",{1.7,2.7}}]

Answer 1

A simply way to obtain the desired plots is to excluding portions of the curves in the discontinuous region. Precisely where t is discontinuous does not seem well defined, so I have taken it to be approximately the almost circular region shown in the first plot in the question. Then,

nodiscontinuity = ContourPlot[{t[x, y] == 2.3, t[x, y] == 2.5, t[x, y] == 2.7}, 
    {x, -1, 1}, {y, -1, 1}, PlotRange -> {{-1, 1}, {-1, 1}}, PlotPoints -> 20,
    MaxRecursion -> 1, 
    RegionFunction -> Function[{x, y, z}, (x - .57)^2 + (y - .01)^2 > .08]]

Superimposing this plot on the result of

ListDensityPlot[data, DataRange -> {{-1, 1}, {-1, 1}}, Evaluate@plot];

Show[%, nodiscontinuity]

then shows that the curves indeed lie just outside the apparent discontinuous region.

(Note that plot (b) in the Question does not lie entirely outside the discontinuous region.)

In answer to the other two questions, the error messages occur because ContourPlot is calling t with symbolic input before calling it with numerical input. And, replacing

data = Table[t[x, y], {y, -1, 1.0, 0.05}, {x, -1, 1.0, 0.05}];

by

data = Table[t[x, y], {x, -1, 1.0, 0.05}, {y, -1, 1.0, 0.05}];

rotates the plot produced by ListDensityPlot, because the second line of code produces the Transpose of the array produced by the first line of code. This can be undone by

ListDensityPlot[data // Transpose, DataRange -> {{-1, 1}, {-1, 1}}, Evaluate@plot]

reproducing the original, not rotated plot.

Addendum

In comments below, Qi Zhong asked whether his creative solution above could be computed more rapidly. This can be accomplished by calling ContourPlot fewer times. However, the 2.1 contour still must be computed separately, because it requires much higher resolution. On my PC, the final plot in the Question require 230 sec, but

Show[ContourPlot[t[x, y], {x, -1, 1}, {y, -1, 1}, PlotPoints -> 20, 
    Contours -> DeleteCases[Range[1.8, 2.6, 0.1], 2.1], ColorFunction -> (White &), 
    Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}, 
    ContourStyle -> (Directive[color[#], Thickness[.005], Opacity[1]] & /@ 
    Range[1.8, 2.6, 0.1]), MaxRecursion -> 1], 
    ContourPlot[t[x, y] == 2.1, {x, -0.33, 0.85}, {y, -0.33, 0.33}, 
    ContourStyle -> color[2.1], PlotPoints -> 20, MaxRecursion -> 1, 
    Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}]]

requires only 106 sec.