First let's save the plot in cp:
cp = ContourPlot3D[
p1 p2^6 (c0 p1 - (-1 + c0) p2)^3 ((-1 + c0) p1 - c0 p2)^2 ==
(-1 + p1) (-1 + p2)^6 (-1 + p1 - (-1 + c0) p2)^3 (1 + (-1 + c0) p1 - c0 p2)^2,
{p1, 0, 0.3}, {p2, 0.5, 1}, {c0, 0, 1},
Lighting -> ({"Directional", White, #} & /@ Tuples[{-1, 1}, 3]),
Mesh -> None, ContourStyle -> {Directive[Red, Opacity[0.5]]},
AxesLabel -> {"\!\(\*SubscriptBox[\(p\), \(1\)]\)",
"\!\(\*SubscriptBox[\(p\), \(2\)]\)", "C"}];
[By the way, you should avoid using capital C, which is used by Mathematica for constant parameters.]
The easy way
Project the coordinates of the polygons in cp onto the p1/p2 plane:
Graphics[
Cases[cp, GraphicsComplex[pts_, g_, rest___] :>
GraphicsComplex[
Most /@ pts, (* the first two coordinates are p1, p2 *)
{Directive[Lighter[Red, 0.5]], Cases[g, _Polygon, Infinity]}
],
Infinity],
Method -> {"TransparentPolygonMesh" -> True}]
I replaced the semi-opacity with a semi-light red. This should be fine if you keep the graphics in the background. Otherwise, it's harder.
Graphics styling.
To get axes on the graph and other features, add Axes -> True and other options:
Graphics[
Cases[cp, GraphicsComplex[pts_, g_, rest___] :>
GraphicsComplex[
Most /@ pts, (* the first two coordinates are p1, p2 *)
{Directive[Lighter[Red, 0.5]], Cases[g, _Polygon, Infinity]}
],
Infinity],
Method -> {"TransparentPolygonMesh" -> True},
Axes -> True]
(See above.) Others to consider are AxesOrigin -> {0, 0}, using Frame -> True instead of Axes, and setting the ratio height/width with AspectRatio -> 1/GoldenRatio.
Edit - The hard way -- Envelopes
First, call the equation eqn and change C to c0:
eqn = (p1 p2^6 (C p1 - (-1 + C) p2)^3 ((-1 + C) p1 - C p2)^2 ==
(-1 + p1) (-1 + p2)^6 (-1 + C p1 - (-1 + C) p2)^3 (1 + (-1 + C) p1 - C p2)^2) /. C -> c0
p1 p2^6 (c0 p1 - (-1 + c0) p2)^3 ((-1 + c0) p1 - c0 p2)^2 ==
(-1 + p1) (-1 + p2)^6 (-1 + c0 p1 - (-1 + c0) p2)^3 (1 + (-1 + c0) p1 - c0 p2)^2
The boundaries of the projection will be either the edges of the 3D plot or the envelope of the surface. The edges of the 3D graph are where c is 0 or 1. There are many solutions which we and we have to extract the ones corresponding to the plot. Luckily the only ones in the plot region are the boundaries and they can be written in the form p2 = func[p1]; hence, the boundaries can picked out by inspecting the value of p2 at some intermediate value of p1, say p1 = 0.15.
sols = Extract[#, Position[p2 /. # /. p1 -> 0.15, p_Real /; 0.5 <= p <= 1]] &@
Solve[eqn && (0 == c0 || c0 == 1), {c0, p2}]
{{c0 -> 1, p2 -> (-1 + p1 + Sqrt[p1 - p1^2])/(-1 + 2 p1)},
{c0 -> 0, p2 -> 1 - p1 ...omitted...}}
Here we have to set up the differential equation for the envelope, which will be the points on the surface where the partial derivative with respect to c0 is zero (envelope == 0); the condition "on the surface" is modeled by the integral curve being orthogonal to the gradient (surface == 0). One of the boundary values of the differential equation is {p1, p2, c0} == {0, 1, 0.} (I've change the value of c0 three times because I've been fooled by the Mathematica output, the numerics being so tricky). The DE is degenerate there, so we have to find a good initial point for NDSolve. As with sols there are many solutions and we have to extract the initial point we want. There is only one envelope line in the plot region, so finding it presents no trouble. The solutions found will be of the form p2 = fp2[p1], c0 = fc0[p1], where fp2, fc0 will be the InterpolatingFunction returned by NDSolve.
fnsrule = Thread[{p2, c0} -> {fp2[p1], fc0[p1]}];
envelope = D[Subtract @@ eqn, {c0}] /. fnsrule;
surface = (D[Subtract @@ eqn, {{p1, p2, c0}}] /. fnsrule) . {1, fp2'[p1], fc0'[p1]};
initsol = NSolve[D[Subtract @@ eqn, {c0}] == 0 && eqn == 0 /. p1 -> 1/100, {p2, c0}];
initpt = First @ Extract[initsol,
Position[{p2, c0} /. initsol, {p_Real, c_Real} /; 0.5 <= p <= 1 && 0 <= c <= 1]]
{p2 -> 0.751159, c0 -> 0.267272}
Expect warning messages. The choice of p1 -> 1/100 was through trial and error. (It's far enough away from p1 == 0 to stay out of serious trouble and give a good solution. Getting much closer tended to call for increasing WorkingPrecision, but increasing WorkingPrecision tended to cause the kernel to crash.) As you can see below, the DE can be solved for almost all the requested domain.
bdy = NDSolve[{envelope == 0, surface == 0,
fp2[0.01] == p2 /. initpt, fc0[0.01] == c0 /. initpt},
{fp2, fc0}, {p1, 0, 0.3}]
{{fp2 -> InterpolatingFunction[{{1.91025*10^-9, 0.3}}, <>],
fc0 -> InterpolatingFunction[{{1.91025*10^-9, 0.3}}, <>]}}
However, the problem still has one more tricky pitfall lurking: The envelope ("green") appears to be the lower boundary of the projection and to be tangent at about p1 == 0.2 to the red line, which corresponds to the boundary of the ContourPlot where c0 == 0. In fact, the envelope passes out of the plot region into a region where c0 < 0, and the red line becomes the lower boundary, which can be seen by inspecting the value of c0 at p1 = 0.3.
boundariesPlot = Plot[Evaluate[Append[p2 /. sols, fp2[p1] /. First@bdy]], {p1, 0, 0.3},
PlotStyle -> {Blue, Red, Darker@Green}]
fc0[0.3] /. First@bdy
-0.488191
So we need to find where the envelope crosses c0 == 0:
tanPt = p1 /. FindRoot[fc0[p1] /. bdy, {p1, 0.2}]
0.196357
Then we can plot the projection of the region: The first solution in sols (where c0 == 1) is the upper boundary. The lower is equal either fp2 or the second solution in sols depending on p1; you can use Piecewise or an If statement.
pl = With[{fp2 = fp2 /. First@bdy},
Plot[
{p2 /. First@sols,
Piecewise[
{{fp2[p1], fp2["Domain"][[1, 1]] <= p1 <= tanPt},
{p2 /. Last@sols, True}}]
},
{p1, 0, 0.3}, Filling -> {1 -> {2}},
FillingStyle -> Directive[Red, Opacity[0.5]], PlotStyle -> None,
PlotRange -> {0.5, 1}]]
You can make a Polygon of the projection from a plot, if you omit PlotStyle -> None. It may be more convenient to combine it with other graphics. Extract the two graphs with Cases and stitch the points together, reversing the direction of one of the lines.
envPts = Join[First[#], Reverse@Last[#]] &@Cases[
With[{fp2 = fp2 /. First@bdy},
Plot[Evaluate[
{p2 /. First@sols,
Piecewise[
{{fp2[p1], fp2["Domain"][[1, 1]] <= p1 <= tanPt},
{p2 /. Last@sols, True}}]
}],
{p1, 0, 0.3}]], Line[pts_] :> pts, Infinity];
envPlot = Graphics[{Opacity[0.5], Red, Polygon[envPts]},
Axes -> True, PlotRange -> {.5, 1}, AspectRatio -> 1/GoldenRatio];
Show[envPlot, boundariesPlot]
Here we can compare the easy and hard ways. The easy way is recolored Blue and so differences will show up as Blue or light Red. Close inspection shows that the easy way has greater granularity typical of surface plots.
Show[proj /. _RGBColor -> Blue, envPlot,
AspectRatio -> 1/GoldenRatio, PlotRange -> {0.5, 1}]
Show[proj /. _RGBColor -> Blue, envPlot,
AspectRatio -> 1, PlotRange -> {{0, 0.02}, {.7, 1}}]
While some steps might be automated, it seems that in any particular case, you're going to have to inspect what's going on and make an intelligent choice. The "easy way" above is fairly straightforward.
Edit 2 - 3rd way - MeshFunctions
We can use MeshFunctions to find the envelope. It can be a good way to solve for a curve on a surface.
This way still has its hangups. The graph gets extremely flat near p1 == 0 and it's hard to compute anything. Still, we can do well enough if we increase the precision and subdivide the plot region as we approach p1 == 0. The following goes down to p1 = 3*10^(-10) an order of magnitude at a time.
Table[
mymesh[i] = With[{prec = 200},
ContourPlot3D[Evaluate @ eqn,
{p1, 3*10^(-i - 1), 3*10^(-i)}, {p2, SetPrecision[0.5, prec], 1}, {c0, 0, 1},
PlotPoints -> 50, MaxRecursion -> 3, Mesh -> {{0}},
BoundaryStyle -> None, ContourStyle -> None,
MeshFunctions -> {Function[{p1, p2, c0},
Evaluate[CubeRoot@D[Subtract @@ eqn, {c0}]]]},
WorkingPrecision -> prec]
], {i, 9}]; // Timing
{4.208199, Null}
The curve is still a ways from reaching p2 == 1, but it doesn't matter for the sake of plotting. We are concerned with p2 as a function of p1, and we don't need p1 much below 10^-4 or 10^-5 -- we can just connect a straight line to the endpoint. Here the segments of the envelope are colored differently:
With[{i1 = 1, i2 = 9},
Graphics3D[
MapThread[Cases[#1,
GraphicsComplex[pts_, g_, rest___] :> GraphicsComplex[pts, {#2, Thick,
Select[Cases[g, _Line, Infinity],
Last @ pts[[ #[[1, 1]] ]] < 0.5 &&
(First@pts[[ #[[1, 1]] ]] > 0.1 || Last@pts[[ #[[1, 1]] ]] > 0.001) &]
}]] &,
{Table[mymesh[i], {i, i1, i2}], Table[Hue[i/9], {i, i1, i2}]}],
PlotRange -> {{0, 0.2}, {0.6, 1}, {0, 0.3}}, BoxRatios -> {1, 1, 1},
Axes -> True,
AxesLabel -> {"\!\(\*SubscriptBox[\(p\), \(1\)]\)",
"\!\(\*SubscriptBox[\(p\), \(2\)]\)", "C"}]
]
Collect the points down to p1 == 3 * 10^-5:
envelopePoints = Append[With[{i1 = 1, i2 = 4},
Flatten[
Map[Cases[#1, GraphicsComplex[pts_, g_, rest___] :>
pts[[
Select[Cases[g, _Line, Infinity],
Last@pts[[#[[1, 1]]]] < 0.5 &&
(First@pts[[#[[1, 1]]]] > 0.1 ||
Last@pts[[#[[1, 1]]]] > 0.001) &][[1, 1]]
]] ] &,
Table[mymesh[i], {i, i1, i2}]],
2]
], {0., 1., 0.}];
Then join them with the other parts of the boundary and compare with the other methods.
envMeshPlot =
Graphics[{Opacity[0.5], Red,
Polygon @ Join[Table[{p1, p2 /. Last@sols}, {p1, 0.3, 0.2, -0.01}],
Most /@ envelopePoints,
Table[{p1, p2 /. First@sols} /. p1 -> p^2, {p, 0, Sqrt[0.3], Sqrt[0.3]/100}]]}]
Show[proj /. _RGBColor -> Blue, envMeshPlot,
AspectRatio -> 1, PlotRange -> {{0, 0.005}, {.7, 1}}, Axes -> False]
Show[envMeshPlot, envPlot /. _RGBColor -> Blue,
AspectRatio -> 1, PlotRange -> {{0, 0.005}, {.7, 1}}, Axes -> False]
The first method is easiest; the drawback is that the overlapping polygons make Opacity give an nonuniform appearance.
RegionPlot[]? $\endgroup$RegionPlot[]works only for flat plane. $\endgroup$