Is there any non-experimental way to obtain CL, CD, and CM as a function of the angle of attack and the control surface deflection angle?
The graphs in the OP show fairly typical results that one would expect from linear analyses (e.g. VLM). 

1. Lift

In the linear range, the lifting surface's lift coefficient ($C_L$) can be expressed as:

$$C_L = a_0 \alpha + a_1 \delta$$

where $a_0=\frac{\partial{C_L}}{\partial{\alpha}}$ is the lift curve slope and $a_1=\frac{\partial{C_L}}{\partial{\delta}}$ is the lift slope per plain-flap deflection, which are constants in the linear range and can be readily seen in your lift plot.

Conclusion: You can simplify your lookup tables to the above equations (all constants) if you are only interested in capturing the linear effect (which is all that VLM can capture anyway).

2. Pitching moment

In the linear range, you can express the pitching moment coefficient ($C_m$), which I assume is computed at quarter chord, as:

$$C_{m} = C_{m_{ac}}(\delta) + C_L\frac{l_{ac}(\delta)}{\overline{c}}$$

where $C_{m_{ac}}$ is the pitching moment at the surface's aerodynamic centre (AC) and is a function of $\delta$; $l_{ac}$ is the distance between the 1/4c and AC and is also a function of $\delta$; $\overline{c}$ is the reference chord length.

The $C_m$ plots you've shown seems to line up very close to the aerodynamic centre of the lifting surface (by the fact that there is very little change in  $C_m$ with respect to AOA). Flap deflection only significantly changes the offset ($C_{m_{ac}}$). The variations in the slopes are because you do not have a full span plain-flap (you have elevons which are partial span), which changes the spanwise lift distribution and the resulting aerodynamic centre.

Conclusion: You can replace your multi-dimensional lookup table to two 1-dim tables for $C_{m_{ac}}$ and $l_{ac}$, and supplant with the above equation. You may even replace these lookup tables with constant linear relationships if your plots hold water.

3. Drag

The drag coefficient ($C_D$) is the most interesting. It seems like there is some boundary layer analysis combined with the VLM, since $C_{D_0}$ changes per flap. If the Reynolds number do not change significantly, you can attempt the following simplication:

$$C_D = C_{D_0}(\delta) + K(\delta)C_L^2$$

where C_{D_0} is the form drag and $K$ is the induced drag factor, both of which are a function of $\delta$.

Conclusion: You may again simplify your multi-dim lookup tables to two 1-dim tables.

4. Lastly...

You definitely will start to see nonlinear effects at 30 deg control deflection in real life, maybe even at 20 deg. This is obviously not captured by VLM.