How to use the splipy.Surface function in Splipy

crv = constructor()
                    elif crv_type == 100:
                        crv = self.splines(1)
                    else:
                        raise IOError('Unsopported trimming curve or malformed input file')
                    two_curves.append(crv)

                # only keep the parametric version (re-generate physical one if we need it)
                one_loop.append(two_curves[0])
                self.trimming_curves.append(two_curves[1])
            all_loops.append(one_loop)

        return TrimmedSurface(surface.bases[0], surface.bases[1], surface.controlpoints, surface.rational, all_loops, raw=True)

    g2_type = [100, 200, 700] # curve, surface, volume identifiers
    classes = [Curve, Surface, Volume]
    g2_generators = {120:line, 130:circle, 140:ellipse,
                     260:cylinder, 292:disc, 270:sphere, 290:torus, 250:plane,
                     210:bounded_surface, 261:surface_of_linear_extrusion} #, 280:cone

    def __init__(self, filename):
        if filename[-3:] != '.g2':
            filename += '.g2'
        self.filename = filename
        self.trimming_curves = []

    def __enter__(self):
        return self

    def write(self, obj):
        if not hasattr(self, 'fstream'):
            self.onlywrite = True

def controlpoints(spline):
    """ Return controlpoints according to nutils ordering """
    n = len(spline)
    dim = spline.dimension
    if isinstance(spline, Curve):
        return np.reshape(spline[:,:]                  , (n, dim), order='F')
    elif isinstance(spline, Surface):
        return np.reshape(spline[:,:,:].swapaxes(0,1)  , (n, dim), order='F')
    elif isinstance(spline, Volume):
        return np.reshape(spline[:,:,:,:].swapaxes(0,2), (n, dim), order='F')
    raise RuntimeError('Non-spline argument detected')

def coons_patch(bottom, right, top, left):
    # coons patch (https://en.wikipedia.org/wiki/Coons_patch)
    top.reverse()
    left.reverse()

    # create linear interpolation between opposing sides
    s1 = edge_curves(bottom, top)
    s2 = edge_curves(left, right)
    s2.swap()
    # create (linear,linear) corner parametrization
    linear = BSplineBasis(2)
    rat = s1.rational  # using control-points from top/bottom, so need to know if these are rational
    if rat:
        bottom = bottom.clone().force_rational() # don't mess with the input curve, make clone
        top.force_rational()                     # this is already a clone
    s3 = Surface(linear, linear, [bottom[0], bottom[-1], top[0], top[-1]], rat)

    # in order to add spline surfaces, they need identical parametrization
    Surface.make_splines_identical(s1, s2)
    Surface.make_splines_identical(s1, s3)
    Surface.make_splines_identical(s2, s3)

    result = s1
    result.controlpoints += s2.controlpoints
    result.controlpoints -= s3.controlpoints
    return result

raise RuntimeError('finitestrain_patch not supported for rational splines')

    # these are given as a oriented loop, so make all run in positive parametric direction
    top.reverse()
    left.reverse()

    # in order to add spline surfaces, they need identical parametrization
    Curve.make_splines_identical(top, bottom)
    Curve.make_splines_identical(left, right)

    # create an initial mesh (correct corners) which we will morph into the right one
    p1 = bottom.order(0)
    p2 = left.order(0)
    p  = max(p1,p2)
    linear = BSplineBasis(2)
    srf = Surface(linear, linear, [bottom[0], bottom[-1], top[0], top[-1]])
    srf.raise_order(p1-2, p2-2)
    for k in bottom.knots(0, True)[p1:-p1]:
        srf.insert_knot(k, 0)
    for k in left.knots(0, True)[p2:-p2]:
        srf.insert_knot(k, 1)

    # create computational mesh
    n1 = len(bottom)
    n2 = len(left)
    dim= left.dimension
    domain, geom = mesh.rectilinear(srf.knots())
    ns = function.Namespace()
    ns.basis = domain.basis('spline', degree(srf), knotmultiplicities=multiplicities(srf)).vector( 2 )
    ns.phi   = domain.basis('spline', degree(srf), knotmultiplicities=multiplicities(srf))
    ns.eye       = np.array([[1,0],[0,1]])
    ns.cp        = controlpoints(srf)

self.scale   = self.width*(1-2*self.margin) / (boundingbox[2]-boundingbox[0])
            self.height  = self.width*geometryRatio + 2*marginPixels
        self.center      = [boundingbox[0], boundingbox[1]]
        self.offset      = [marginPixels, marginPixels]

        # create xml root tag
        self.xmlRoot = etree.Element('svg',  {'xmlns':'http://www.w3.org/2000/svg',
                                              'version':'1.1',
                                              'width':str(self.width),
                                              'height':str(self.height)})

        # populate tree with all curves and surfaces in entities
        for entry in self.all_objects:
            if isinstance(entry, Curve):
                self.write_curve(self.xmlRoot, entry)
            elif isinstance(entry, Surface):
                self.write_surface(entry)

        # if no objects are stored, then we've most likely only called read()
        if len(self.all_objects) > 0:
            rough_string = etree.tostring(self.xmlRoot) # entire xml-file on one line
            reparsed     = minidom.parseString(rough_string)
            result       = reparsed.toprettyxml(indent="  ") # adds newline and inline
            f = open(self.filename, 'w')
            f.write(result)

m = len(circle_seg) # number of control points of the sweep
    cp = np.zeros((m * n, 4))

    # loop around the circle and set control points by the traditional 9-point
    # circle curve with weights 1/sqrt(2), only here C0-periodic, so 8 points
    dt = 0
    t  = 0
    for i in range(m):
        x,y,w = circle_seg[i]
        dt  = atan2(y,x) - t
        t  += dt
        curve.rotate(dt)
        cp[i * n:(i + 1) * n, :]  = curve[:]
        cp[i * n:(i + 1) * n, 2] *= w
        cp[i * n:(i + 1) * n, 3] *= w
    result = Surface(curve.bases[0], circle_seg.bases[0], cp, True)
    # rotate it back again
    result.rotate(normal_phi,   [0,1,0])
    result.rotate(normal_theta, [0,0,1])
    return result

"""
    type = kwargs.get('type', 'coons')
    if len(curves) == 1: # probably gives input as a list-like single variable
        curves = curves[0]
    if len(curves) == 2:
        crv1 = curves[0].clone()
        crv2 = curves[1].clone()
        Curve.make_splines_identical(crv1, crv2)
        (n, d) = crv1.controlpoints.shape  # d = dimension + rational

        controlpoints = np.zeros((2 * n, d))
        controlpoints[:n, :] = crv1.controlpoints
        controlpoints[n:, :] = crv2.controlpoints
        linear = BSplineBasis(2)

        return Surface(crv1.bases[0], linear, controlpoints, crv1.rational)
    elif len(curves) == 4:
        # reorganize input curves so they form a directed loop around surface
        rtol = state.controlpoint_relative_tolerance
        atol = state.controlpoint_absolute_tolerance
        mycurves = [c.clone() for c in curves] # wrap into list and clone all since we're changing them
        dim = np.max([c.dimension for c in mycurves])
        rat = np.any([c.rational  for c in mycurves])
        for i in range(4):
            for j in range(i+1,4):
                Curve.make_splines_compatible(mycurves[i], mycurves[j])

        if not (np.allclose(mycurves[0][-1], mycurves[1][0], rtol=rtol, atol=atol) and
                np.allclose(mycurves[1][-1], mycurves[2][0], rtol=rtol, atol=atol) and
                np.allclose(mycurves[2][-1], mycurves[3][0], rtol=rtol, atol=atol) and
                np.allclose(mycurves[3][-1], mycurves[0][0], rtol=rtol, atol=atol)):
            reorder = [mycurves[0]]

:return: The utah teapot
    :rtype: List of Surface
    """
    path = join(dirname(realpath(__file__)), 'templates', 'teapot.bpt')
    with open(path) as f:
        results = []
        numb_patches = int(f.readline())
        for i in range(numb_patches):
            p = np.fromstring(f.readline(), dtype=np.uint8, count=2, sep=' ')
            basis1 = BSplineBasis(p[0]+1)
            basis2 = BSplineBasis(p[1]+1)

            ncp = basis1.num_functions() * basis2.num_functions()
            cp  = [np.fromstring(f.readline(), dtype=np.float, count=3, sep=' ') for j in range(ncp)]
            results.append(Surface(basis1, basis2, cp))

    return results

[ -4./3*(2*sr3-1),               0, 4./3*(1-2*sr3),   4*(5-sr3)/3], # row 2
              [-sr2/3*(7-2*sr3),               0,       -5*sr6/3, sr2*(sr3+6)/3],
              [              0 ,               0,    4*(sr3-5)/3, 4*(5*sr3-1)/9],
              [ sr2/3*(7-2*sr3),               0,       -5*sr6/3, sr2*(sr3+6)/3],
              [  4./3*(2*sr3-1),               0, 4./3*(1-2*sr3),   4*(5-sr3)/3],
              [    -sr2*(4-sr3),             sr2,    sr2*(sr3-4), sr2*(3*sr3-2)], # row 3
              [    -(3*sr3-2)/2,    -(2-3*sr3)/2,     -(sr3+6)/2,     (sr3+6)/2],
              [              0 ,-sr2*(2*sr3-7)/3,       -5*sr6/3, sr2*(sr3+6)/3],
              [     (3*sr3-2)/2,    -(2-3*sr3)/2,     -(sr3+6)/2,     (sr3+6)/2],
              [     sr2*(4-sr3),             sr2,    sr2*(sr3-4), sr2*(3*sr3-2)],
              [      -4*(sr3-1),      -4*(1-sr3),      4*(1-sr3),   4*(3-sr3)  ], # row 4
              [           -sr2 ,    -sr2*(sr3-4),    sr2*(sr3-4), sr2*(3*sr3-2)],
              [              0 , -4./3*(1-2*sr3), 4./3*(1-2*sr3),4./3*(5-sr3)  ],
              [            sr2 ,    -sr2*(sr3-4),    sr2*(sr3-4), sr2*(3*sr3-2)],
              [       4*(sr3-1),      -4*(1-sr3),      4*(1-sr3),   4*(3-sr3)  ]]
        wmin = Surface(b,b,cp, rational=True)
        wmax = wmin.clone().mirror([0,0,1])
        vmax = wmin.clone().rotate(pi/2, [1,0,0])
        vmin = vmax.clone().mirror([0,1,0])
        umax = vmin.clone().rotate(pi/2, [0,0,1])
        umin = umax.clone().mirror([1,0,0])
        # ideally I would like to call edge_surfaces() now, but that function
        # does not work with rational surfaces, so we'll just manually try
        # and add some inner controlpoints
        cp   = np.zeros((5,5,5,4))
        cp[ :, :, 0,:] = wmin[:,:,:]
        cp[ :, :,-1,:] = wmax[:,:,:]
        cp[ :, 0, :,:] = vmin[:,:,:]
        cp[ :,-1, :,:] = vmax[:,:,:]
        cp[ 0, :, :,:] = umin[:,:,:]
        cp[-1, :, :,:] = umax[:,:,:]
        inner = np.linspace(-.5,.5, 3)