'Elliptic solver for bending energy minimization.
Uses a dedicated \'toy\' sparse Jacobi PCG solver.'

def compute_dz(self):

    dz_init = _DOF_estimator_geom.compute_dz(self)
    Uf0 = np.ravel(dz_init)
    reference_element = _ReducedHCT_Element()
    J = CubicTriInterpolator._get_jacobian(self._tris_pts)
    eccs = self._eccs
    triangles = self._triangles
    Uc = self.z[self._triangles]
    (Kff_rows, Kff_cols, Kff_vals, Ff) = reference_element.get_Kff_and_Ff(J, eccs, triangles, Uc)
    tol = 1e-10
    n_dof = Ff.shape[0]
    Kff_coo = _Sparse_Matrix_coo(Kff_vals, Kff_rows, Kff_cols, shape=(n_dof, n_dof))
    Kff_coo.compress_csc()
    (Uf, err) = _cg(A=Kff_coo, b=Ff, x0=Uf0, tol=tol)
    err0 = np.linalg.norm((Kff_coo.dot(Uf0) - Ff))
    if (err0 < err):
        warnings.warn(u'In    TriCubicInterpolator    initialization,    PCG    sparse    solver    did    not    converge    after    1000    iterations.    `geom`    approximation    is    used    instead    of    `min_E`')
        Uf = Uf0
    dz = np.empty([self._pts.shape[0], 2], dtype=np.float64)
    dz[:, 0] = Uf[::2]
    dz[:, 1] = Uf[1::2]
    return dz
