diff --git a/manimlib/mobject/mobject.py b/manimlib/mobject/mobject.py
index 2df95afb..48059363 100644
--- a/manimlib/mobject/mobject.py
+++ b/manimlib/mobject/mobject.py
@@ -22,7 +22,7 @@ from manimlib.utils.simple_functions import get_parameters
 from manimlib.utils.space_ops import angle_of_vector
 from manimlib.utils.space_ops import get_norm
 from manimlib.utils.space_ops import rotation_matrix_transpose
-from manimlib.utils.shaders import ShaderWrapper
+from manimlib.shader_wrapper import ShaderWrapper
 
 
 # TODO: Explain array_attrs
diff --git a/manimlib/mobject/types/surface.py b/manimlib/mobject/types/surface.py
index 9c8e4a1f..94455ced 100644
--- a/manimlib/mobject/types/surface.py
+++ b/manimlib/mobject/types/surface.py
@@ -44,6 +44,7 @@ class ParametricSurface(Mobject):
         self.uv_func = uv_func
         self.compute_triangle_indices()
         super().__init__(**kwargs)
+        self.sort_faces_back_to_front()
 
     def init_points(self):
         dim = self.dim
@@ -124,7 +125,7 @@ class ParametricSurface(Mobject):
         self.rgbas = interpolate(mobject1.rgbas, mobject2.rgbas, alpha)
         return self
 
-    def pointwise_become_partial(self, smobject, a, b):
+    def pointwise_become_partial(self, smobject, a, b, axis=1):
         assert(isinstance(smobject, ParametricSurface))
         self.points[:] = smobject.points[:]
         if a <= 0 and b >= 1:
@@ -132,25 +133,38 @@ class ParametricSurface(Mobject):
 
         nu, nv = smobject.resolution
         self.points[:] = np.vstack([
-            self.get_partial_points_array(arr, a, b, (nu, nv, 3))
+            self.get_partial_points_array(arr, a, b, (nu, nv, 3), axis=axis)
             for arr in self.get_surface_points_and_nudged_points()
         ])
         return self
 
-    def get_partial_points_array(self, points, a, b, resolution):
+    def get_partial_points_array(self, points, a, b, resolution, axis):
         nu, nv = resolution[:2]
         points = points.reshape(resolution)
-        lower_index, lower_residue = integer_interpolate(0, nv - 1, a)
-        upper_index, upper_residue = integer_interpolate(0, nv - 1, b)
-        tuples = [
-            (points[:, :lower_index], lower_index, lower_residue),
-            (points[:, upper_index:], upper_index, upper_residue),
-        ]
-        for to_change, index, residue in tuples:
-            col = interpolate(points[:, index], points[:, index + 1], residue)
-            to_change[:] = np.tile(col, to_change.shape[1]).reshape(to_change.shape)
+        max_index = resolution[axis] - 1
+        lower_index, lower_residue = integer_interpolate(0, max_index, a)
+        upper_index, upper_residue = integer_interpolate(0, max_index, b)
+        if axis == 0:
+            points[:lower_index] = interpolate(points[lower_index], points[lower_index + 1], lower_residue)
+            points[upper_index:] = interpolate(points[upper_index], points[upper_index + 1], upper_residue)
+        else:
+            tuples = [
+                (points[:, :lower_index], lower_index, lower_residue),
+                (points[:, upper_index:], upper_index, upper_residue),
+            ]
+            for to_change, index, residue in tuples:
+                col = interpolate(points[:, index], points[:, index + 1], residue)
+                to_change[:] = col.reshape((nu, 1, *resolution[2:]))
         return points.reshape((nu * nv, *resolution[2:]))
 
+    def sort_faces_back_to_front(self, vect=OUT):
+        tri_is = self.triangle_indices
+        indices = list(range(len(tri_is) // 3))
+        indices.sort(key=lambda i: np.dot(self.points[tri_is[3 * i]], vect))
+        for k in range(3):
+            tri_is[k::3] = tri_is[k::3][indices]
+        return self
+
     # For shaders
     def get_shader_data(self):
         s_points, du_points, dv_points = self.get_surface_points_and_nudged_points()
@@ -241,13 +255,13 @@ class TexturedSurface(ParametricSurface):
                 sm.set_opacity(opacity, family)
         return self
 
-    def pointwise_become_partial(self, tsmobject, a, b):
-        super().pointwise_become_partial(tsmobject, a, b)
+    def pointwise_become_partial(self, tsmobject, a, b, axis=1):
+        super().pointwise_become_partial(tsmobject, a, b, axis)
         self.im_coords[:] = tsmobject.im_coords
         if a <= 0 and b >= 1:
             return self
         nu, nv = tsmobject.resolution
-        self.im_coords[:] = self.get_partial_points_array(self.im_coords, a, b, (nu, nv, 2))
+        self.im_coords[:] = self.get_partial_points_array(self.im_coords, a, b, (nu, nv, 2), axis)
         return self
 
     def get_shader_uniforms(self):
diff --git a/manimlib/utils/shaders.py b/manimlib/shader_wrapper.py
similarity index 100%
rename from manimlib/utils/shaders.py
rename to manimlib/shader_wrapper.py
diff --git a/manimlib/shaders/quadratic_bezier_fill_vert.glsl b/manimlib/shaders/quadratic_bezier_fill_vert.glsl
index 469c1e7f..dc426d63 100644
--- a/manimlib/shaders/quadratic_bezier_fill_vert.glsl
+++ b/manimlib/shaders/quadratic_bezier_fill_vert.glsl
@@ -6,7 +6,6 @@ uniform float is_fixed_in_frame;
 in vec3 point;
 in vec3 unit_normal;
 in vec4 color;
-// in float fill_all;  // Either 0 or 1
 in float vert_index;
 
 out vec3 bp;  // Bezier control point