_clough_tocher_2d_single direction between centroids

[tludwigsmr]

@pv

Hi Pauli,

The way the cross-edge derivatives are enforced in the scipy implementation of the Clough-Tocher interpolant seems very reasonable (although I do not understand yet all details). However, there is a small issue:

When there is no neighbour triangle, according to the comments, the (non-existing) neighbour centroid V_4' is assumed on the line (e_12 + e_13)/2 (e.g. at the mid-edge), and the coefficients g1, g2, g3 should be calculated such that V_4'-V_4 = const_[(V_4-V_2) + g_i_(V_3 - V_2)].

It looks like in this case, the values for g1,g2,g3 should be rather -1/2 than -2/3.

This can be checked by taking an equilateral triangle with the following coordinates, values, and slopes:

coor = [(0, 0), (1, 0), (0.5, math.cos(math.radians(30)))]
values = [1, 0, 0]
slopes = [(0, 0), (0, 0), (0, 0)]

With a value of -2/3, the interpolant is not symmetric w.r.t. (e_12 + e_13)/2, as shown below:

N.B.: I have not verified this with the scipy code, but with a stand-alone Python script where I compare different HCT- and other interpolants, including an (almost) 1:1 copy of _clough_tocher_2d_single.

Best regards,
Thomas

[pv]

Thanks, I tried to recheck the result now and also got -1/2. I don't remember the details of how this works anymore either, but it seems there's a calculation error.

The implementation probably would have been more clear if it was done in the FEM way by affinely mapping prototype elements, but I didn't at the time think about the problem too deeply...