How do I transpose an ellipse function to stretch the ellipse into curved space?

Suppose your parabola is given by $$y=a x^2+b x + c$$ Its derivative is $y'=2 a x + b$. A vector perpendicular to your parabola at the point $(x,y)$ is thus $$(-2 a x - b,1)$$ Let us normalize this vector $$\vec{n}=\frac{(-2 a x - b,1)}{\sqrt{(2 a x + b)^2+1}}$$ The curves at distance $d$ from your parabola (green in your figure) are thus $$(x,a x^2+b x + c)+d\vec{n}$$ The above formula transforms cartesian-like coordinates $(x,d)$ into the transformed space $$T(x,d)=(x,a x^2+b x + c)+d\vec{n}$$ All we need to do is feed the ellipse expression to this transformation $$T(r_1\cos(\theta)+r_1,r_2\sin(\theta)+r_2)$$ this gives a curve with the desired characteristics.

Example:

This curve touches the straight lines at their midpoints. It also touces the top and bottom curve once, but not at their exact midpoint. If you want a curve that touches the bottom curves at its precise midpoint, you need to parametrize your parabola by its arc length, which is not possible in terms of elementary functions.

As others have noted, this is really just a deformation or change-of-coordinates problem. Our ellipse is defined in a rectangular cartesian coordinate system, and we want to map it to another system where coordinates are distance along the parabola and normal distance away from the parabola.

These kinds of deformations are common in high-end CAD systems. However, as @Wouter pointed out, the equation of the deformed ellipse is the composition of the original ellipse equation and the change-of-variables transformation, so it's an ugly thing that no CAD system can represent exactly (as far as I know). So, the CAD system will typically give you a spline approximation of the true deformed shape. This is perfectly OK, of course, as long as the approximation is more accurate than your manufacturing tolerances.

The picture below shows the result I got in my favorite CAD system.I didn't use the same numbers as you, but that doesn't matter. The point is to see what can be achieved with standard built-in functions in typical CAD systems.

More info in these answers.

The map $$s\mapsto\left\{\eqalign{x_0(s)&:={\rm arsinh}\, s \cr y_0(s)&:=1-\sqrt{1+s^2}\cr}\right.\qquad(-\infty<s<\infty)$$ maps the $s$-axis isometrically onto the catenary $$\gamma:\qquad y=1-\cosh x=-{1\over2}x^2-{1\over24}x^4+?x^6\ ,$$ which is in the intended range a very good approximation to the parabola $y=-{1\over2}x^2$. One computes $$\dot x_0(s)={1\over\sqrt{1+s^2}},\qquad\dot y_0(s)=-{s\over\sqrt{1+s^2}}\ ,$$ confirming that $\dot x_0^2(s)+\dot y_0^2(s)\equiv1$. It follows that the map $$(s,t)\mapsto\left\{\eqalign{x(s,t)&:={\rm arsinh}\, s+{t s\over\sqrt{1+s^2}} \cr y(s,t)&:=1-\sqrt{1+s^2}+{t\over\sqrt{1+s^2}}\cr}\right.\qquad(-\infty<s<\infty, \ -h<t<h)\tag{1}$$ maps the horizontal strip $|t|<h$, $h>0$ sufficiently small, of the $(s,t)$-plane bijectively onto a curved strip along the catenary $\gamma$ in such a way that (i) the soul $t=0$ is mapped isometrically onto $\gamma,$ and (ii) vertical segments $t\mapsto (s_0,t)$ $(-h<t<h)$ are mapped isometrically onto vertical segments orthogonal to $\gamma$ at $\bigl(x_0(s_0),y_0(s_0)\bigr)\in\gamma$. In this way horizontal lines $s\mapsto(s, t_0)$ $(-\infty<s<\infty)$ are mapped onto parallel curves of $\gamma$ at distance $|t_0|$ from $\gamma$.

If you now are given a curve $$\theta\mapsto\bigl(s(\theta),t(\theta)\bigr)\qquad(0\leq\theta\leq\Theta)\tag{2}$$ in the $(s,t)$-plane, e.g., an ellipse, then you immediately obtain a parametric representation of its distorted image in the $(x,y)$-plane by composing the representation $(2)$ with $(1)$.