I have a function that is expensive to evaluate whose inputs are n-dimensional (n is the order of a dozen or two).  I need the output of this function at each node and each time step for a PDE simulation.  I know that the portion of the n-d space sampled is small since n-1 of the inputs are species mole fractions.  I know the outputs are continuous in all dimensions.  I would like to have the solver adaptively sample the n-d space as needed based on the inputs it is being provided.  If the new inputs are close to/contained by old inputs, a simple linear interpolation would be sufficient.  If the new inputs are significantly different from old inputs, I want to add points to the table.  The inputs/outputs are relatively compact, so RAM is not an issue.  I am just wondering if there are data structures/algorithms that will make the lookup and interpolation efficient while allowing for the quick addition of new points.
Adaptive Table Lookup for Expensive Function Evaluation
interpolation