How to use the matscipy.fracture_mechanics.crack.CubicCrystalCrack function in matscipy

def test_consistency_of_deformation_gradient_and_stress(self):
        for C11, C12, C44, surface_energy, k1 in self.materials:
            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
            k = crack.k1g(surface_energy)*k1
            for i in range(10):
                x = np.random.uniform(-10, 10)
                y = np.random.uniform(-10, 10)
                F = crack.deformation_gradient(x, y, 0, 0, k)
                eps = (F+F.T)/2-np.eye(2)
                sig_x, sig_y, sig_xy = crack.stresses(x, y, 0, 0, k)
                eps_xx = crack.crack.a11*sig_x + \
                         crack.crack.a12*sig_y + \
                         crack.crack.a16*sig_xy
                eps_yy = crack.crack.a12*sig_x + \
                         crack.crack.a22*sig_y + \
                         crack.crack.a26*sig_xy
                # Factor 1/2 because elastic constants and matrix product are
                # expressed in Voigt notation.
                eps_xy = (crack.crack.a16*sig_x + \

refcell = None
            reftip_x = None
            reftip_y = None
            #[41, 39, 1], 
            for i, n in enumerate([[21, 19, 1], [11, 9, 1], [6, 5, 1]]):
                #print(n)
                cryst = diamond(el, a0, n, crack_surface, crack_front)
                set_groups(cryst, n, skin_x, skin_y)
                cryst.set_pbc(True)
                cryst.set_calculator(calc)
                FIRE(UnitCellFilter(cryst), logfile=None).run(fmax=1e-6)
                cryst.set_cell(cryst.cell.diagonal(), scale_atoms=True)

                ase.io.write('cryst_{}.xyz'.format(i), cryst, format='extxyz')

                crk = crack.CubicCrystalCrack(crack_surface,
                                              crack_front,
                                              Crot=C/units.GPa)
                k1g = crk.k1g(surface_energy)

                tip_x = cryst.cell.diagonal()[0]/2
                tip_y = cryst.cell.diagonal()[1]/2

                a = cryst.copy()
                a.set_pbc([False, False, True])

                k1 = 1.0
                ux, uy = crk.displacements(cryst.positions[:,0], cryst.positions[:,1],
                                           tip_x, tip_y, k1*k1g)

                a.positions[:, 0] += ux
                a.positions[:, 1] += uy

def test_elastostatics(self):
        eps = 1e-3
        for C11, C12, C44, surface_energy, k1 in self.materials:
            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
            k = crack.k1g(surface_energy)*k1
            for i in range(10):
                x = i+1
                y = i+1
                # Finite difference approximation of the stress divergence
                sxx0x, syy0x, sxy0x = crack.stresses(x-eps, y, 0, 0, k)
                sxx0y, syy0y, sxy0y = crack.stresses(x, y-eps, 0, 0, k)
                sxx1x, syy1x, sxy1x = crack.stresses(x+eps, y, 0, 0, k)
                sxx1y, syy1y, sxy1y = crack.stresses(x, y+eps, 0, 0, k)
                divsx = (sxx1x-sxx0x)/(2*eps) + (sxy1y-sxy0y)/(2*eps)
                divsy = (sxy1x-sxy0x)/(2*eps) + (syy1y-syy0y)/(2*eps)
                # Check that divergence of stress is zero (elastostatic
                # equilibrium)
                self.assertAlmostEqual(divsx, 0.0, places=3)
                self.assertAlmostEqual(divsy, 0.0, places=3)

return

        nx = 128
        for calc, a0, C11, C12, C44, surface_energy, bulk_coordination in [
            ( atomistica.DoubleHarmonic(k1=1.0, r1=1.0, k2=1.0,
                                        r2=math.sqrt(2), cutoff=1.6),
              clusters.sc('He', 1.0, [nx,nx,1], [1,0,0], [0,1,0]),
              3, 1, 1, 0.05, 6 ),
#            ( atomistica.Harmonic(k=1.0, r0=1.0, cutoff=1.3, shift=False),
#              clusters.fcc('He', math.sqrt(2.0), [nx/2,nx/2,1], [1,0,0],
#                           [0,1,0]),
#              math.sqrt(2), 1.0/math.sqrt(2), 1.0/math.sqrt(2), 0.05, 12)
            ]:
            print '{} atoms.'.format(len(a0))

            crack = CubicCrystalCrack(C11, C12, C44, [1,0,0], [0,1,0])

            x, y, z = a0.positions.T
            r2 = min(np.max(x)-np.min(x), np.max(y)-np.min(y))/4
            r1 = r2/2

            a = a0.copy()
            a.center(vacuum=20.0, axis=0)
            a.center(vacuum=20.0, axis=1)
            ref = a.copy()
            r0 = ref.positions

            a.set_calculator(calc)
            print 'epot = {}'.format(a.get_potential_energy())

            sx, sy, sz = a.cell.diagonal()
            tip_x = sx/2

def test_J_integral(self):
        if quadrature is None:
            print('No scipy.integrate.quadrature. Skipping J-integral test.')
            return

        for C11, C12, C44, surface_energy, k1 in self.materials:
            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
            k = crack.k1g(surface_energy)*k1

            def polar_path(theta, r=1, x0=0, y0=0):
                nx = np.cos(theta)
                ny = np.sin(theta)
                n = np.transpose([nx, ny])
                return r*n-np.array([x0, y0]), n, r

            def elliptic_path(theta, r=1, x0=0, y0=0):
                rx, ry = r
                x = rx*np.cos(theta)
                y = ry*np.sin(theta)
                nx = ry*np.cos(theta)
                ny = rx*np.sin(theta)
                ln = np.sqrt(nx**2+ny**2)
                nx /= ln

def test_consistency_of_deformation_gradient_and_displacement(self):
        eps = 1e-3
        for C11, C12, C44, surface_energy, k1 in self.materials:
            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
            k = crack.k1g(surface_energy)*k1
            for i in range(10):
                x = i+1
                y = i+1
                F = crack.deformation_gradient(x, y, 0, 0, k)
                # Finite difference approximation of deformation gradient
                u, v = crack.displacements(x, y, 0, 0, k)
                ux0, vx0 = crack.displacements(x-eps, y, 0, 0, k)
                uy0, vy0 = crack.displacements(x, y-eps, 0, 0, k)
                ux1, vx1 = crack.displacements(x+eps, y, 0, 0, k)
                uy1, vy1 = crack.displacements(x, y+eps, 0, 0, k)
                du_dx = (ux1-ux0)/(2*eps)
                du_dy = (uy1-uy0)/(2*eps)
                dv_dx = (vx1-vx0)/(2*eps)
                dv_dy = (vy1-vy0)/(2*eps)
                du_dx += np.ones_like(du_dx)

for nx in [ 8, 16, 32, 64 ]:
            for calc, a, C11, C12, C44, surface_energy, bulk_coordination in [
                #( atomistica.DoubleHarmonic(k1=1.0, r1=1.0, k2=1.0,
                #                            r2=math.sqrt(2), cutoff=1.6),
                #  clusters.sc('He', 1.0, [nx,nx,1], [1,0,0], [0,1,0]),
                #  3, 1, 1, 0.05, 6 ),
                (
                  #atomistica.Harmonic(k=1.0, r0=1.0, cutoff=1.3, shift=True),
                  IdealBrittleSolid(k=1.0, a=1.0, rc=1.3),
                  clusters.fcc('He', math.sqrt(2.0), [nx,nx,1], [1,0,0],
                               [0,1,0]),
                  math.sqrt(2), 1.0/math.sqrt(2), 1.0/math.sqrt(2), 0.05, 12)
                ]:
                clusters.set_groups(a, (nx, nx, 1), 0.5, 0.5)
                crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)

                a.center(vacuum=20.0, axis=0)
                a.center(vacuum=20.0, axis=1)
                a.set_calculator(calc)

                sx, sy, sz = a.cell.diagonal()
                tip_x = sx/2
                tip_y = sy/2

                k1g = crack.k1g(surface_energy)
                r0 = a.positions.copy()

                u, v = crack.displacements(a.positions[:,0], a.positions[:,1],
                                           tip_x, tip_y, k1g)
                a.positions[:,0] += u
                a.positions[:,1] += v