I have a finite element solver that I implemented in MATLAB. I am calculating a specific potential function that is constant in each tetrahedral element. The question is, I want to be able to interpolate this potential on a rectangular grid on a slice. How can I achieve this? I can take the middle points of each element and then utilize a point-by-point interpolation. However, I also want to utilize the information that the potential is constant inside the element itself.

I want this interpolation to calculate the directional derivatives of the potential. It would be even better to be able to calculate directional derivatives of this potential. I cannot do these inside elements, obviously, since the potential is constant in each element. If you could suggest a method to calculate directional partial derivatives, it would be even better.

  • 2
    $\begingroup$ I'm not quite clear about the details of the question because you don't distinguish between the original solution and the interpolant when you say "potential". Can you add formulas in which you give each of these names? $\endgroup$ Feb 15, 2021 at 14:53

1 Answer 1


I did something similar for reporting stress in a linear elasticity finite element program, which approximates stress as a constant quantity across a linear tetrahedral element. It is often desirable to plot stress as a continuous quantity, because a smoothly varying quantity is more believable. I figured this was a good compromise. This approach should be valid for whatever quantity you wish to interpolate.

In short, I interpolated the stress at each of the vertices, while still honoring the elemental average values. For each element e, approximate the value at each node a corresponding to that element using your basis functions $$ \sigma_e = \sum_a N_a\sigma_a d\Gamma $$ Integrate over the element. $$ \int_e\sigma_e d\Gamma = \int_e\ \sum_a N_a\sigma_a d\Gamma $$ Equivalently. $$ \sigma_eV_e = \sum_a (\int_e N_a d\Gamma)\sigma_a $$ Therefore, you end up with an underdetermined system of e equations and v variables, where e is the number of elements and v is the number of vertices. An underdetermined system results in a least squares approximation that honors each of your equations. Technically, the nodes at the boundary are a linear extrapolation of that data, but the average quantities are still honored.

Here is a plot of the original "real" data (stress is a symmetric 3x3 matrix -- this is pressure which is a scalar reduction of that quantity) that is constant at each element.

Pressure -- Elemental Values

Here is a plot of the interpolated "nodal" data that varies linearly between vertices.

Pressure -- Nodal Values

EDIT: I only just saw that you wanted to take directional derivatives. However, I think you can do this as you now have a linear approximation of the quantity inside the element, which is certainly differentiable.

  • $\begingroup$ Thank you very much for your great answer and please forgive me for not being able to read it earlier. I will definitely try this method and share the results with you. Cheers! $\endgroup$
    – strahd
    Feb 20, 2021 at 12:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.