'self.df is computed as weighted average of _triangles sharing a common
node. On each triangle itri f is first assumed linear (= ~f), which
allows to compute d~f[itri]
Then the following approximation of df nodal values is then proposed:
f[ipt] = SUM ( w[itri] x d~f[itri] , for itri sharing apex ipt)
The weighted coeff. w[itri] are proportional to the angle of the
triangle itri at apex ipt'

def compute_dz(self):

    el_geom_w = self.compute_geom_weights()
    el_geom_grad = self.compute_geom_grads()
    w_node_sum = np.bincount(np.ravel(self._triangles), weights=np.ravel(el_geom_w))
    dfx_el_w = np.empty_like(el_geom_w)
    dfy_el_w = np.empty_like(el_geom_w)
    for iapex in range(3):
        dfx_el_w[:, iapex] = (el_geom_w[:, iapex] * el_geom_grad[:, 0])
        dfy_el_w[:, iapex] = (el_geom_w[:, iapex] * el_geom_grad[:, 1])
    dfx_node_sum = np.bincount(np.ravel(self._triangles), weights=np.ravel(dfx_el_w))
    dfy_node_sum = np.bincount(np.ravel(self._triangles), weights=np.ravel(dfy_el_w))
    dfx_estim = (dfx_node_sum / w_node_sum)
    dfy_estim = (dfy_node_sum / w_node_sum)
    return np.vstack([dfx_estim, dfy_estim]).T
