Computation of shape functions
Shape functions are approximations of linear partial differential operators in given points.
Derivation
Let $\mathcal{L}$ be a linear partial differential operator. We wish to approximate $(\mathcal{L}u)(\vec{p})$ using only values of $u$ in support of $p$. To do so, we construct a MLS approximation (see Moving Least Squares (MLS)) and compute the coefficients $\b{\alpha}$, such that $\hat{u}(\vec{p}) = \b{b}(\vec{p})^\T \b{\alpha}$. The coefficients are computed by solving a system $WB \b{\alpha} = W\b{u}$ using the pseudoinverse \[ \b{\alpha} = (WB)^+W\b{u}, \] where $A^+$ denotes the Moore-Penrose pseudoinverse that can be calculated using SVD decomposition. Writing down the approximation function $\hat{u}$ we get \[ \hat u (\vec{p}) = \b{b}(\vec{p})^\T \b{\alpha} = \b{b}(\vec{p})^\T (WB)^+W\b{u} = \b{\chi}_{\vec{p}} \b{u}. \]
We have defined $\b{\chi}$ to be \[ \b{\chi}(\vec{p}) = \b{b}(\vec{p})^\T (WB)^+W. \] Vector $\b{\chi}$ is a row vector, also called a shape function. The name comes from being able to take all the information about shape of the domain and choice of approximation and store it in a single row vector, being able to approximate a function value from given support values $\b{u}$ with a single dot product. For any values $\b{u}$, value $\b{\chi}(p) \b{u}$ gives us the approximation $\hat{u}(\vec{p})$ of $u$ in point $\vec{p}$. Mathematically speaking, $\b{\chi}(\vec{p})$ is a functional, $\b{\chi}(\vec{p})\colon \R^n \to \R$, mapping $n$-tuples of known function values to their approximations in point $\vec{p}$.
The same approach works for any linear operator $\mathcal L$ applied to $u$. We approximate \[ (\mathcal{L}u)(\vec{p}) \approx (\mathcal{L} \hat{u})(p) = (\mathcal{L}\b{b})(\vec{p})^\T \b{\alpha} = (\mathcal{L}\b{b})(\vec{p})^\T (WB)^+W\b{u} = \b{\chi}_{\mathcal{L},\vec{p}} \b{u}. \]
For example, take a $1$-dimensional case for approximation of derivatives with weight equal to $1$ and $n=m=3$, with equally spaced support values at distances $h$. We wish to approximate $u''$ in the middle support point, just by making a weighted sum of the values, something like the finite difference \[ u'' \approx \frac{u_1 - 2u_2 + u_3}{h^2}. \] This is exactly the same formula as we would have come to by computing $\b{\chi}$, except that our approach is a lot more general. But one should think about $\b{\chi}$ as one would about the finite difference scheme, it is a rule, telling us how to compute the derivative. \[ u''(s_2) \approx \underbrace{\begin{bmatrix} \frac{1}{h^2} & \frac{-2}{h^2} & \frac{1}{h^2} \end{bmatrix}}_{\b{\chi}} \begin{bmatrix}u_1 \\ u_2 \\ u_3 \end{bmatrix} \]
The fact that $\b{\chi}$ is independent of the function values $\b{u}$ but depend only on domain geometry means that we can just compute the shape functions $\b{\chi}$ for points of interest and then approximate any linear operator of any function, given its values, very fast, using only a single dot product.
Numerical calculation of the shape functions
The expression \[ \b{\chi}_{\mathcal{L}, \vec{p}} = (\mathcal{L}\b{b})(\vec{p})^\T (WB)^+W \] can be evaluated directly, but this is not the most optimal approach. Numerically cheaper and more stable way is to translate the problem of inverting the matrix to solving a linear system of equations.
Invertible $B$ case: If $B$ is invertible, then $\b{\chi}_{\mathcal{L}, \vec{p}} = (\mathcal{L}\b{b})(\vec{p})^\T B^{-1}$, transposing the equation then multiplying it from the left by $B$, $\b{\chi}_{\mathcal{L}, \vec{p}}$ can be thought as a solution of a system $B^\T\b{\chi}_{\mathcal{L}, \vec{p}}^\T = (\mathcal{L}\b{b})$, which can be solved using LU or Cholesky decomposition for example.
General case: For a system written as $Ax = b$, where $A$ is a $n\times m$ matrix, $x$ is a vector of length $m$ and $b$ a vector of length $n$, a generalized solution is defined as $x$ that minimizes $\|A x - b\|_2^2$. If more $x$ attain the minimal value, $x$ with the minimal $\|x\|$ is chosen. Note that this generalizes the solution a general system ($A$ is invertible) and over-determined system ($n > m$ and $A$ has full rank). Such an $x$ can be computed using the pseudoinverse $x = A^{+} b$.
In our case, let us denote a part of the solution containing the pseudoinverse by $\tilde{\b{\chi}}_{\mathcal{L}, \vec{p}}$. \[ \b{\chi}_{\mathcal{L}, \vec{p}}(\vec{p}) = \underbrace{(\mathcal{L}\b{b})(\vec{p})^\T (WB)^+}_{\tilde{\b{\chi}}_{\mathcal{L}, \vec{p}}} W \] We have an expression $\tilde{\b{\chi}}_{\mathcal{L}, \vec{p}} = (\mathcal{L}\b{b})(\vec{p})^\T (WB)^+$ which after transposition takes the form $\tilde{\b{\chi}}_{\mathcal{L}, \vec{p}}^\T = ((WB)^\T)^+(\mathcal{L}\b{b})(\vec{p})$, the same as $x = A^+b$ above. Therefore, $\tilde{\b{\chi}}_{\mathcal{L}, \vec{p}}^\T$ is the solution of an (indeterminate) system $(WB)^\T \tilde{\b{\chi}}_{\mathcal{L}, \vec{p}}^\T = (\mathcal{L}\b{b})(\vec{p})$. After solving that, we can get the shape function $\b{\chi}_{\mathcal{L}, \vec{p}} = \tilde{\b{\chi}}_{\mathcal{L}, \vec{p}} W$ by multiplying by matrix $W$.
The system before can be solved using any decomposition of matrix $(WB)^\T = B^\T W$, not neceserally the SVD decompostion, but depending on our knowledge of the problem, we can use Cholesky ($B^\T W$ is positive definite), $LDL^\T$ if it is symmetric, $LU$ for a general square matrix, $QR$ for full rank overdetermined system and SVD for a general system. If more shapes need to be calculated using the same matrix $B^\T W$ and only different right hand sides, it can be done efficiently by storing the decomposition of $B^\T W$.
