'Parameters
*alpha* is a (N x 3 x 1) array (array of column-matrices) of
barycentric coordinates
*J* is a (N x 2 x 2) array of jacobian matrices (jacobian matrix at
triangle first apex)
*ecc* is a (N x 3 x 1) array (array of column-matrices) of triangle
eccentricities
*dofs* is a (N x 1 x 9) arrays (arrays of row-matrices) of computed
degrees of freedom.
Returns
Returns the values of interpolated function 2nd-derivatives
[d2z/dx2, d2z/dy2, d2z/dxdy] in global coordinates at locations alpha,
as a column-matrices of shape (N x 3 x 1).'

def get_function_hessians(self, alpha, J, ecc, dofs):

    d2sdksi2 = self.get_d2Sidksij2(alpha, ecc)
    d2fdksi2 = _prod_vectorized(dofs, d2sdksi2)
    H_rot = self.get_Hrot_from_J(J)
    d2fdx2 = _prod_vectorized(d2fdksi2, H_rot)
    return _transpose_vectorized(d2fdx2)
