Which mathematical subfield or formalism deals with local but non-pointwise function approximations? Main Question
I have a class of high-dimensional functions which obey certain constraints (of value or partial derivative) on specific submanifolds of their domain.
I am quite confident that there is a more specific class of functions that matches these functions on these submanifolds and approximates them in the vicinity.
I want to prove this or find out what additional assumptions I need to make.
To this end, I am looking for a mathematical subfield or formalism that deals with such questions.
Example
A simple statement along the lines of what I wish to prove is:

If $f(0,y) = f(x,0) = 0$ and $f$ is sufficiently benign, there exist functions $g$, $h$, and $r$ such that:
  
  
*
  
*$f(x,y) = g(x) h(y) + r(x,y)$,
  
*$g(0) = h(0) = 0$,
  
*$r(x,0) = r(0,y) = 0$,
  
*$∂_1 r(0,y) = ∂_2 r(x,0) = 0$.
  

In this case, the answer is to choose
$g(x) = \frac{∂_2 f(x,0)}{∂_1 ∂_2 f(0,0)}$,
$h(x) = ∂_1 f(0,y)$, and
$r(x,y) = f(x,y) - g(x) h(y)$.
What I considered so far
Taylor series allow me to perform pointwise function approximations.
If done sufficiently general, this may allow me to find the shape of the desired approximation, but the process is rather tedious and suggests that there is a more general variant of it.
Approximation theory seems to focus on numerically finding best approximations on the entire domain without paying special attention to certain regions – which is not what I want.
 A: You can prove such statements employing the Mather division theorem.
For example, you can use the following statement derived from the theorem:

Let $f$, $f_1$ be smooth functions satisfying the conditions of the example and let additionally $∂_x ∂_y f_1(0,0)\ne0$. Then there exist a neighbourhood of the origin $U$ and a smooth function $q$ s.t. $f(x,y)=q(x,y)f_1(x,y)$ in $U$.

To see that, we'll make an orthogonal change of coordinates 
$u=(x+y)/2$,$v=(x−y)/2$, so $xy=u^2−v^2$, and apply the Mather division theorem.
Overloading notation we have: $f_1(0,0)=\partial_uf_1(0,0)$. To get $k=2$ in the theorem one needs $\partial_u^2f_1(0,0)\ne0$. But this derivative is a linear combination of $\partial_x^2f_1(0,0)$, $\partial_y^2f_1(0,0)$ and 
$\partial_x\partial_yf_1(0,0)$. Since the first two are equal to zero, the condition $\partial_u^2f_1(0,0)\ne0$ is equivalent to  $∂_x ∂_y f_1(0,0)\ne0$.
By the Mather division theorem:
$$f(u,v)=q(u,v)f_1(u,v)+r_0(v)+r_1(v)u.$$
From 
$$\begin{alignat*}{1}
f(v,v)=r_0(v)+r_1(v)v=0,\\
f(-v,v)=r_0(v)-r_1(v)v=0,\end{alignat*}$$
it follows that $r_0≡r_1≡0$ and $f(x,y)=q(x,y)f_1(x,y)$. 
Taking $f_1(x,y)=xy$, one has $f(x,y)=q(x,y)xy$ and the solution for your example can be rewritten as $g(x)h(y)=\frac{xyq(x,0)q(0,y)}{q(0,0)}$.
