For some physics simulations on surfaces, it can be useful (for efficiency, simplicity, etc.) to first simulate an object in 2D as a sector before wrapping into 3D. For example, if you are interested in modelling a sector of a cloth draped over a conical surface. The 3D object can be thought of as being generated by rolling the 2D sector around its origin.

In these situations, it’s necessary to compute the pushforward for the 2D vector field to 3D space, an operation from differential geometry. This post breaks down the process for doing so, and for pulling the vectors back, which might be useful for post processing analysis.

General Remark: While I focus on the example of rolling a 2D sector onto a cone, everything here is actually a special case of mapping from a 2D parameter domain \((r,\phi)\) into a 3D surface \((x(r,\phi),\,y(r,\phi),\,z(r,\phi))\). In more general settings, we have a surface \(S \subset \mathbb{R}^3\) parameterized by \(F(u,v) = (x(u,v), y(u,v), z(u,v))\). The pushforward of a 2D vector \(\langle{v^u, v^v}\rangle{}\) is given by \(v^{u}\frac{\partial F}{\partial u} + v^{v}\frac{\partial F}{\partial v}\). The same approach—deriving the Jacobian \(\frac{\partial(x,y,z)}{\partial(r,\phi)}\) and then using it to push vectors forward—applies to any smooth parametric surface, whether it’s a cylinder, a sphere (in patches), or something more complicated.

We start with a velocity vector defined in a local 2D Cartesian plane and need to map it into a 3D Cartesian system. This means we need to transform its positional coordinates \((x,y)\) to \((x,y,z)\) and its component vectors \((v_x, v_y)\) to \((v_x,v_y,v_z)\).

However, instead of going directly from \((v_x,v_y)\) to \((v_x,v_y,v_z)\), we will first switch to polar coordinates \((r,\phi)\) and components \((v_r, v_{\phi})\) to take advantage of natural cylindrical symmetry, and then embed into 3D.

Step 1: From \((x,y)\) to \((r,\phi)\) and \((v_x,v_y)\) to \((v_r, v_{\phi})\)

The vector lives at position \((x,y)\) on the sector. First, let’s get the polar coordinates of the vector’s position (NOT to be confused with the radial and tangential components of the vector).

import numpy as np
# Generate random example x,y coords
num_samples = 10 
xy = np.random.uniform(-100, 100, (num_samples, 2)) 
x = xy[:, 0]
y = xy[:, 1]
r = np.sqrt(x**2 +y**2) # or np.linalg.norm(xy)
phi = np.arctan2(y, x) # arctan(y,x)

Although we didn’t enforce it here, let’s pretend these \((x,y)\) points are certainly on the sector of interest.

Next, to map a vector \((v_x, v_y)\) to its polar form \((v_r, v_{\phi})\), we need to express it in terms of its radial and tangential components. \(\begin{cases} v_r = v_xcos(\phi)+ v_ysin(\phi), \\ v_{\phi} = -v_xsin(\phi) + v_ycos(\phi) \end{cases}\)

\(v_r\) tells us how much of the vector is increasing in the direction \(r\), and \(v_{\phi}\) tells us how much of the vector is increasing in the direction perpendicular to \(r\), i.e. tangent to the circular motion.

# Random example vector components
vxy = np.random.uniform(-10,10, (num_samples, 2)) 
vx = vxy[:, 0]
vy = vxy[:, 1]
v_r  = vx * np.cos(phi) + vy * np.sin(phi)
v_phi = -vx * np.sin(phi) + vy * np.cos(phi)

Side Note on Coordinate vs. Orthonormal Vectors: Here, \((v_r, v_{\phi})\) are the coordinate-basis components in polar coordinates. In some applications (e.g for certain PDEs on surfaces), you might want an orthonormal basis \(\hat{r}, \hat{\phi}\), Then \(\hat{\phi} = (\frac{1}{r})\frac{\partial}{\partial \phi}\), and the formula for \(v_{\phi}\) would involve different factors of \(r\). For this post, we stick to coordinate-basis compnents – just keep in mind if you need unit vectors in the future.

Now that we’ve got the vector’s original coordinates and components in polar space, we’re ready to embed the vector in 3D space.

Step 2: From \((r,\phi)\) to \((x,y,z)\)

First, note that since we are rolling a periodic sector of a circle into a cone, we need to obtain symmetry factor \(k\), which describes the number of times the cone needs to be wrapped in order to follow the periodicity. \(k\) is \(2\pi\) (the ‘full angle’ in radians) divided by the opening angle of the sector \(\alpha\), so \(k = \frac{2\pi}{\alpha}\). This means if, for example, \(\alpha = \frac{\pi}{2}\), then \(k = 4\), and the intuition behind that is that the sector would be wrapped four times around the cone. Note that the derivation of the Jacobian below is only valid for \(\alpha < 2\pi\).

alpha = np.pi/(np.random.uniform(1, 4)) # let's say alpha is somewhere between pi/4 and pi
k = 2.0 * np.pi/alpha 
theta = phi * k # New angular coordinate on the cone

Tip / Boundary Caution: At \(r = 0\) then we are at the cone tip. The Jacobian can degenerate there. Practically, that just means that you have to be careful pushing forward or back at exactly at the tip; the geometry has a singularity. You will see when we introduce an epsilon factor to deal with this later. In other applications, it can be fine to exclude \(r=0\) from the domain or handle it with special boundary conditions depending on the real world analog.

We need the Jacobian of \((x,y,z)\) with respect to \((r,\phi)\) to see how infinitesimal displacements \((\partial{r}, \partial{\phi})\) in 2D become displacements in \((dx, dy, dz)\). The Jacobian acts as a change-of-basis matrix that allows us to correctly express how vectors in the 2D domain correspond to vectors in 3D. Or in other words, to account for how distances and angles change under the rolling process, which introduces non-linear distortions.

Using the transformation equations:

We compute the partial derivatives:

1. Partial derivatives with respect to \(r\): \(\frac{\partial x}{\partial r} = \frac{1}{k} \cos(\theta), \quad \frac{\partial y}{\partial r} = \frac{1}{k} \sin(\theta), \quad \frac{\partial z}{\partial r} = \frac{1}{k} \sqrt{k^2 - 1}.\)

2. Partial derivatives with respect to \(\phi\): \(\frac{\partial x}{\partial \phi} = -\frac{r}{k} \sin(\theta) \cdot k = -r \sin(\theta),\) \(\frac{\partial y}{\partial \phi} = \frac{r}{k} \cos(\theta) \cdot k = r \cos(\theta),\) \(\frac{\partial z}{\partial \phi} = 0.\)

Thus, the Jacobian matrix is:

df_dr = np.stack(
    (
        np.cos(theta) / k,
        np.sin(theta) / k,
        np.full_like(theta, np.sqrt(k**2 - 1) / k)
    ),
    axis = -1
)

df_dphi = np.stack(
    (
        -r * np.sin(theta),
        np.cos(theta),
        np.zeros_like(theta)
    ),
    axis = -1
)

v_xyz = v_r[..., np.newaxis] * df_dr + v_phi[..., np.newaxis] * df_dphi

import numpy as np 
def pushforward(
    xy: np.ndarray,
    vxy: np.ndarray,
    alpha: float,
) -> np.ndarray: 
    k = 2.0 * np.pi/alpha 

    r = np.sqrt(xy[...,0]**2 +xy[...,1]**2) # or np.linalg.norm(xy)
    phi = np.arctan2(xy[...,1], xy[...,0]) # arctan(y,x)
    theta = phi * k  # New angular coordinate on the cone
    
    # polar components of vector
    v_r  = vxy[..., 0] * np.cos(phi) + vxy[..., 1] * np.sin(phi)
    v_phi = -vxy[..., 0] * np.sin(phi) + vxy[..., 1] * np.cos(phi)


    df_dr = np.stack(
        (
            np.cos(theta) / k,
            np.sin(theta) / k,
            np.full_like(theta, np.sqrt(k**2 - 1) / k)
        ),
        axis = -1
    )

    df_dphi = np.stack(
        (
            -r * np.sin(theta),
            np.cos(theta),
            np.zeros_like(theta)
        ),
        axis = -1
    )

    return v_r[..., np.newaxis] * df_dr + v_phi[..., np.newaxis] * df_dphi

Pullback

First we must recover \(r\) and \(\phi\) from the cone point. There is a trap here, however, as \(\phi\) cannot be directly recovered due to the ambiguity introduced by \(\theta = \operatorname{atan2}(y,x)\). The function [\(\operatorname{atan2}\)] returns values in the range \((-\pi, \pi)\), while actual angle \(\theta = k\phi\) is periodic naturally wraps modulo \(2\pi\) on the base of the cone. Because of this, using \(\operatorname{atan2}\) alone would lose information about the true value of \(\phi\), making an accurate pullback impossible.

Global vs Local: This is a prime example of how parameterization can cause trouble if the image “wraps around.” Locally, you can invert \(\theta \to \phi\), but globally, you might need to keep track of extra offset or stored data about which “turn” the unwrapped sector was on.

Why the Transpose?

When we do the pullback, conceputally we are computing: \((v_r, v_{\phi}) = \langle{v_{xyz}, \frac{\partial{F}}{\partial{r}}}\rangle{}, \frac{\partial{F}}{\partial{\phi}}\rangle{}.\) From a linear algebra perspective, that is precisely the transpose of the Jacobian matrix. Hence “pullback = transpose of pushforward” (in coordinate-basis form when the coordinate-basis vectors are orthogonal).

Before we pullback, we need a function that compute \((x,y,z)\) on the cone. In order to do so, all we need are the original \((x,y)\) coordinates and the opening angle \(\alpha\). We borrow from the math earlier that showed this mapping.

def map_xy_to_xyz(
    xy: np.ndarray,
    alpha: float,
) -> np.ndarray:
    k = 2.0 * np.pi/alpha 
    r = np.sqrt(xy[...,0]**2 + xy[...,1]**2)
    phi = np.arctan2(xy[...,1],xy[...,0])
    theta = phi * k 
    # scale down radius by symmetry factor to ensure correct perimeter length 
    x = (r/k) * np.cos(theta)
    y = (r/k) * np.sin(theta)
    z = (r/k) * np.sqrt(k**2 - 1)
    return np.stack([x,y,z], axis=-1)

Great, now we can assume we have the x,y,z coordinates on the cone. We are ready for the pullback.

def pullback(
    xyz: np.ndarray,
    v_xyz: np.ndarray,
    original_phi: np.ndarray,
    alpha: float,
) -> np.ndarray: 

    k = 2.0 * np.pi / alpha 
    theta = original_phi * k 

    # Recover r from the cone coordinates
    cone_x = xyz[..., 0]
    cone_y = xyz[..., 1]
    r = k * np.sqrt(cone_x**2 + cone_y**2)


    df_dr_t = np.stack(
        [
            np.cos(theta)/k,
            np.sin(theta)/k,
            np.full_like(theta, np.sqrt(k**2 - 1)/k),
        ],
        axis = -1,
    )

    df_dphi_t = np.stack(
        [
            -r * np.sin(theta),
            r * np.cos(theta),
            np.zeros_like(theta)
        ],
        axis = -1
    )

    # Recover the polar components of the 2d vector
    epsilon = 1e-8 # include epsilon to avoid divide by 0 at (x = 0, y = 0)
    v_r = (v_xyz * df_dr_t).sum(axis=-1)
    v_phi = (v_xyz * df_dphi_t).sum(axis=-1) / (r**2 + epsilon) 

    # Convert polar components back to Cartesian 
    v_x = np.cos(original_phi) * v_r - np.sin(original_phi) * v_phi
    v_y = np.sin(original_phi) * v_r + np.cos(original_phi) * v_phi

    return np.stack([v_x, v_y], axis = -1) 

Coordinate vs. Orthonormal Again: Depending on how you define \((v_r, v_{\phi})\), you might do “divide by \(r\)” in that last step (rather than \(r^2\)). Our choice is consistent with using the raw coordinate basis for polar coordinates – namely, \(\frac{\partial}{\partial{r}}\) and \(\frac{\partial}{\partial{\phi}}\). In that coordinate basis, the vector \(\frac{\partial}{\partial{\phi}}\) has length \(r\). Consequently, if a vector \(v\) in the plane is written as \(v = v_r\frac{\partial}{\partial{r}} + v_{\phi}\frac{\partial}{\partial{\phi}}\), then \(v_{\phi}\) is the coefficient of a basis vector whose norm is \(r\).

If, instead, you prefer an orthonormal basis \({\hat{r}, \hat{\phi}}\), you would define \(\hat{r} = \frac{\partial}{\partial{r}} / \left\|\frac{\partial}{\partial{r}}\right\|\), \(\hat{\phi} = \frac{\partial}{\partial{\phi}} / \left\|\frac{\partial}{\partial{\phi}}\right\|\). And in standard polar coordinates \(\left\|\frac{\partial}{\partial{\phi}}\right\| = r\). That means \(\hat{\phi} = (\frac{1}{r})\frac{\partial}{\partial{\phi}}.\)

So an orthonormal decomposition would look like \(v = (v_r^{(orth.)})\hat{r} + (v_{\phi}^{(orth.)})\hat{\phi}\), which translates back to the coordinate basis via:

Hence, if you had originally stored \((v_r^{(orth.)}, v_{\phi}^{(orth.)})\) in an orthonormal basis, you would need an extra factor of \(r\) somewhere when converting to (or from) the raw coordinate basis. In short, the difference in dividing by \(r^2\) vs. \(r\) vs. no factor is all about which basis you are using for your 2D vectors:

• Raw coordinate basis: \((\frac{\partial}{\partial{r}}, \frac{\partial}{\partial{\phi}})\).
• Orthonormal basis: \({\hat{r}, \hat{\phi}}\), with \(\hat{\phi} = \frac{1}{r}\frac{\partial}{\partial{\phi}}\).

Our cone example and its code stick to the simpler coordinate basis convention, which is why you see (\(v_{\phi}\))-related operations divided by \(r^2\) or \(r\) in a few places. If you switch to an orthonormal basis, you’ll need to adjust those factors accordingly.

Putting it all together

num_samples = 10 
xy = np.random.uniform(-100, 100, (num_samples, 2)) 
vxy = np.random.uniform(-10,10, (num_samples, 2)) 
alpha = np.pi/(np.random.uniform(1, 4)) # let's say alpha is somewhere between pi/4 and pi

# Push forward
v_xyz = pushforward(xy, vxy, alpha)
#W
# Remember we need to store the original phi
original_phi = np.arctan2(xy[..., 1], xy[..., 0])

xyz = map_xy_to_xyz(xy, alpha)

# Pull back
vxy = pullback(
    xyz, 
    v_xyz, 
    original_phi, 
    alpha
)

Geometric Intuition

• Going via \((r,\theta)\) is useful if the 3D surface naturally uses polar-like coordinates (e.g., cones, cylinders, paraboloids).
• The Jacobian \(\frac{\partial (x,y,z)}{\partial (r,\theta)}\) encodes how motions in \((r,\theta)\) “push forward” into 3D space.
• Rolling/unrolling: The factor \(k = \frac{2\pi}{\alpha}\) ensures that when you go around the cone once, you go around \(2\pi\) in \(\theta\). This matches the sector’s arc length \(\alpha \times r\).
• Tip singularities: At \(r = 0\), the cone’s apex. Watch for degenerate Jacobians.

Beyond Cones: This same approach generalizes to any surface parameterization \((u,v) \to (x(u,v), y(u,v), z(u,v))\). If your surface is a simple patch (like a rectangle), you can still compute partial derivatives w.r.t. \((u,v)\) and push forward vectors or pull them back. The only difference is how you define \((x(u,v), y(u,v), z(u,v))\).

Plots

Pending… (this stuff can lead to some nice plots…)