The following example shows panelization with twisted surfaces. A twisted surface is created by rotating some of control points around other control points.
The sample input file used in the example is below.

surface4.3dm

add_library('igeo')

size(480,360,IG.GL)

IG.open("surface4.3dm")
surfaces = IG.surfaces()

for surf in surfaces :
    unum = 20
    vnum = 20
    uinc = 1.0/unum
    vinc = 1.0/vnum
    for i in range(unum) : 
        for j in range(vnum) : 
            if (i+j)%2== 0 : 
        
                pts = []
                pts.append([])
                pts.append([])
                pts[0].append(surf.pt( i*uinc, j*vinc ))
                pts[0].append(surf.pt( i*uinc, (j+2.0/3)*vinc ))
                pts[0].append(surf.pt( i*uinc, (j+4.0/3)*vinc ))
                pts[0].append(surf.pt( i*uinc, (j+2)*vinc ))
                pts[1].append(surf.pt( (i+1)*uinc, j*vinc ))
                pts[1].append(surf.pt( (i+1)*uinc, (j+2.0/3)*vinc ))
                pts[1].append(surf.pt( (i+1)*uinc, (j+4.0/3)*vinc ))
                pts[1].append(surf.pt( (i+1)*uinc, (j+2)*vinc ))
                
                rotAxis = pts[1][3].dif(pts[1][2])
                
                # rotating around the center pts[1][2] and the axis rotAxis
                pts[0][2].rot(pts[1][2],rotAxis,-PI/2)
                # rotating around the center pts[1][3] and the axis rotAxis
                pts[0][3].rot(pts[1][3],rotAxis,-PI/2)
                
                ISurface(pts, 1, 3)
    surf.del()

Then the rotation angle is randomized.

add_library('igeo')

size(480,360,IG.GL)

IG.open("surface4.3dm")
surfaces = IG.surfaces()

for surf in surfaces : 
    unum = 20
    vnum = 20
    uinc = 1.0/unum
    vinc = 1.0/vnum
    for i in range(unum) : 
        maxAngle = IRand.get(0, PI*2/3)
        for j in range(vnum) : 
            if (i+j)%2== 0 : 
                
                pts = []
                pts.append([])
                pts.append([])
                pts[0].append(surf.pt( i*uinc, j*vinc ))
                pts[0].append(surf.pt( i*uinc, (j+2.0/3)*vinc ))
                pts[0].append(surf.pt( i*uinc, (j+4.0/3)*vinc ))
                pts[0].append(surf.pt( i*uinc, (j+2)*vinc ))
                pts[1].append(surf.pt( (i+1)*uinc, j*vinc ))
                pts[1].append(surf.pt( (i+1)*uinc, (j+2.0/3)*vinc ))
                pts[1].append(surf.pt( (i+1)*uinc, (j+4.0/3)*vinc ))
                pts[1].append(surf.pt( (i+1)*uinc, (j+2)*vinc ))
                
                rotAxis = pts[1][3].dif(pts[1][2])
                
                # rotating around the center pts[1][2] and the axis rotAxis
                pts[0][2].rot(pts[1][2],rotAxis,-maxAngle*j*vinc)
                # rotating around the center pts[1][3] and the axis rotAxis
                pts[0][3].rot(pts[1][3],rotAxis,-maxAngle*j*vinc)
                
                ISurface(pts, 1, 3).clr(0,maxAngle/PI,0)
    surf.del()