# Paint a region with different color in ContourPlot3D

surface[x_, y_, z_] := Piecewise[{{(x – 1/2)^2 + y^2 + z^2 –
1/16, (x – 1/2)^2 + y^2 + z^2 – 1/16 <= 0 && z >= 0}}, z];
ContourPlot3D[{surface[x, y, z] == 0},
{x, -1, 1}, {y, -1, 1}, {z, -0.1, 0.3},
PlotRange -> {{-1, 1}, {-1, 1}, {-0.5, 1.5}},
ImageSize -> 600,
BoundaryStyle -> None,
Mesh -> None,
ViewPoint -> {2, -6, 2}
]

The plot is like below. Then I want to paint the region constrained by

320x^2-304x+256y^2+57>=0 && 15x^2-16x+19y^2+3<=0 with Black color, and the rest area is Orange. ================= ================= 3 Answers 3 ================= This works without making 2 plots and combining them with Show, and it also doesn't require using excessive PlotPoints. I'm not sure why MeshShading works better in this case than ColorFunction, but it does. Read Michael E2's excellent explanation here. The key points are Set the value of Mesh to a single number, 0 Set the MeshFunctions to the two polynomials which you want to be less than and greater than 0, the number you set already. Set the values for MeshShading, giving two values for each MeshFunction, either Automatic or Black, for regions above and below the value of Mesh surface[x_, y_, z_] := Piecewise[{{(x - 1/2)^2 + y^2 + z^2 - 1/16, (x - 1/2)^2 + y^2 + z^2 - 1/16 <= 0 && z >= 0}}, z];
ContourPlot3D[{surface[x, y, z] == 0},
{x, -1, 1}, {y, -1, 1}, {z, -0.1, 0.3},
PlotRange -> {{-1, 1}, {-1, 1}, {-0.5, 1.5}},
ImageSize -> 600,
BoundaryStyle -> None,
Mesh -> {{0}},
ViewPoint -> {2, -6, 2},
MeshFunctions -> {Function[{x, y, z},
320 x^2 – 304 x + 256 y^2 + 57],
Function[{x, y, z}, 15 x^2 – 16 x + 19 y^2 + 3]},
MeshShading -> {{Automatic, Black}, {Black, Automatic}}]

This doesn’t take any longer than the plotting command in the OP.

ColorFunction, I believe, uses VertexColors to apply a color to a vertex of the polygonal mesh that makes the surface of the plot. The vertex colors are interpolated over the interior of a given polygon. The result tends to be blotchy, unless the polygons are very small, and to be computationally somewhat more expensive (in the FE/graphics rendering). Neither ColorFunction nor MeshFunctions factor into mesh refinement, AFAIK, so boundaries in each method are sometimes not smooth.
– Michael E2
May 4 at 10:39

cp1 = ContourPlot3D[{surface[x, y, z] == 0}, {x, -1, 1}, {y, -1,
1}, {z, -0.1, 0.3}, PlotRange -> {{-1, 1}, {-1, 1}, {-0.5, 1.5}},
ImageSize -> 600,
ContourStyle -> Directive[Orange, Opacity[0.8], Specularity[White, 30]],
BoundaryStyle -> None, Mesh -> None, ViewPoint -> {2, -6, 2},
RegionFunction -> Function[{x, y, z}, 320 x^2 – 304 x + 256 y^2 + 57 <=0|| 15 x^2 - 16 x + 19 y^2 + 3>= 0]];

cp2 = ContourPlot3D[{surface[x, y, z] == 0}, {x, -1, 1}, {y, -1,
1}, {z, -0.1, 0.3}, PlotRange -> {{-1, 1}, {-1, 1}, {-0.5, 1.5}},
ImageSize -> 600, ContourStyle -> Black,
RegionFunction -> Function[{x, y, z}, 320 x^2 – 304 x + 256 y^2 + 57 >= 0 &&
15 x^2 – 16 x + 19 y^2 + 3 <= 0], BoundaryStyle -> None, Mesh -> None, ViewPoint -> {2, -6, 2}];

Show[cp1, cp2]

Can it be written in one ContourPlot3D?
– Qi Zhong
May 3 at 20:23

@QiZhong, maybe using ColorFunction. I will post an update if i can make it work.
– kglr
May 3 at 20:25

Thank you. This method is good enough. But there are two problems. First, two ContourPlot3D means the computer needs to do the calculation twice, which needs almost twice time. Second, the intersection region shows part orange, part black when rotating the figure, because there are two layers there. That’s why I want to use only one ContouPlot3D.
– Qi Zhong
May 3 at 20:31

@QiZhong, for the second issue, use RegionFunction -> Function[{x, y, z, f}, 320 x^2 – 304 x + 256 y^2 + 57 <= 0 || 15 x^2 - 16 x + 19 y^2 + 3 >= 0], in cp1.
– kglr
May 3 at 20:40

@kglr I failed to get anything useful from ColorFunction, have you succeeded that way?
– BlacKow
May 3 at 20:52

This can be done with ColorFunction and ColorFunctionScaling.

With

region = 320 x^2 – 304 x + 256 y^2 + 57 >= 0 && 15 x^2 – 16 x + 19 y^2 + 3 <= 0 Then ContourPlot3D[{surface[x, y, z] == 0}, {x, -1, 1}, {y, -1, 1}, {z, -0.1, 0.3}, PlotRange -> {{-1, 1}, {-1, 1}, {-0.5, 1.5}},
ImageSize -> 600, BoundaryStyle -> None, Mesh -> None,
ViewPoint -> {2, -6, 2},
ColorFunctionScaling -> False,
ColorFunction -> Function[{x, y, z, f}, If[region, Black, Orange]]]

You can increase the PlotPoints for a sharper boundary.

Hope this helps.

Yeah, it can work. I need to increase PlotPoints to 40 and MaxRecursion->1. It takes ~5minutes. Is it possible to use plotpoints->20 in the x<0 region and plotpoints0>40 in x>0 region?
– Qi Zhong
May 4 at 1:28