The following example shows panelization with tubed surfaces. A tube surface below is created with 3 control points in u direction and the surface is closed in u direction.
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) : 
            pt1 = surf.pt(i*uinc, j*vinc)
            pt2 = surf.pt((i+1)*uinc, j*vinc)
            pt3 = surf.pt((i+1)*uinc, (j+1)*vinc)
            pt4 = surf.pt(i*uinc, (j+1)*vinc)
            
            depth = 2.0
            pt1d = surf.pt(i*uinc, j*vinc, depth)
            pt2d = surf.pt((i+1)*uinc, j*vinc, depth)
            pt3d = surf.pt((i+1)*uinc, (j+1)*vinc, depth)
            
            cpts1 = []
            cpts1.append([])
            cpts1.append([])
            cpts1.append([])
            cpts1[0].append(pt1)
            cpts1[0].append(pt1d)
            cpts1[1].append(pt2)
            cpts1[1].append(pt2d)
            cpts1[2].append(pt3)
            cpts1[2].append(pt3d)
            
            # u degree 2 (curve), v degree 1 (straight)
            # true for closing surface in u direction
            ISurface(cpts1, 2, 1, True, False)
    surf.del()

In the following code, one tube surface is divided into 4 surfaces connecting the curved profile and the rectangular profile.

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) :
            pt1 = surf.pt(i*uinc, j*vinc)
            pt2 = surf.pt((i+1)*uinc, j*vinc)
            pt3 = surf.pt((i+1)*uinc, (j+1)*vinc)
            pt4 = surf.pt(i*uinc, (j+1)*vinc)
            
            depth = -2.0
            pt1d = surf.pt(i*uinc, j*vinc, depth)
            pt2d = surf.pt((i+1)*uinc, j*vinc, depth)
            pt3d = surf.pt((i+1)*uinc, (j+1)*vinc, depth)
            
            pt1m = pt1d.mid(pt2d)
            pt2m = pt2d.mid(pt3d)
            pt3m = pt3d.mid(pt1d)
            
            # side surface 1
            cpts2 = []
            cpts2.append([])
            cpts2.append([])
            cpts2.append([])
            cpts2[0].append(pt1) 
            cpts2[0].append(pt1m) 
            cpts2[1].append(pt1.mid(pt4))
            cpts2[1].append(pt1d)
            cpts2[2].append(pt4)
            cpts2[2].append(pt3m)
            ISurface(cpts2, 2, 1)
            
            # side surface 2
            cpts3 = []
            cpts3.append([])
            cpts3.append([])
            cpts3.append([])
            cpts3[0].append(pt3)
            cpts3[0].append(pt2m)
            cpts3[1].append(pt3.mid(pt2))
            cpts3[1].append(pt2d)
            cpts3[2].append(pt2)
            cpts3[2].append(pt1m)
            ISurface(cpts3, 2, 1)
            
            # side surface 3
            cpts4 = []
            cpts4.append([])
            cpts4.append([])
            cpts4.append([])
            cpts4[0].append(pt4)
            cpts4[0].append(pt3m)
            cpts4[1].append(pt4.mid(pt3))
            cpts4[1].append(pt3d)
            cpts4[2].append(pt3)
            cpts4[2].append(pt2m)
            ISurface(cpts4, 2, 1)
            
            # side filling triangle
            ISurface(pt1m, pt2, pt1)
    surf.del()

Then the depths of tubular surfaces are 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) :
        for j in range(vnum) :
            pt1 = surf.pt(i*uinc, j*vinc)
            pt2 = surf.pt((i+1)*uinc, j*vinc)
            pt3 = surf.pt((i+1)*uinc, (j+1)*vinc)
            pt4 = surf.pt(i*uinc, (j+1)*vinc)
            
            depth = IRand.get(-5,5)
            dratio = (depth+5)/10
            
            pt1d = surf.pt(i*uinc, j*vinc, depth)
            pt2d = surf.pt((i+1)*uinc, j*vinc, depth)
            pt3d = surf.pt((i+1)*uinc, (j+1)*vinc, depth)
            
            pt1m = pt1d.mid(pt2d)
            pt2m = pt2d.mid(pt3d)
            pt3m = pt3d.mid(pt1d)
            
            cpts2 = []
            cpts2.append([])
            cpts2.append([])
            cpts2.append([])
            cpts2[0].append(pt1)
            cpts2[0].append(pt1m)
            cpts2[1].append(pt1.mid(pt4))
            cpts2[1].append(pt1d)
            cpts2[2].append(pt4)
            cpts2[2].append(pt3m)
            ISurface(cpts2, 2, 1).clr(dratio/3+0.1,0.1,0.1)

            cpts3 = []
            cpts3.append([])
            cpts3.append([])
            cpts3.append([])
            cpts3[0].append(pt3)
            cpts3[0].append(pt2m)
            cpts3[1].append(pt3.mid(pt2))
            cpts3[1].append(pt2d)
            cpts3[2].append(pt2)
            cpts3[2].append(pt1m)
            ISurface(cpts3, 2, 1).clr(dratio/3+0.1,0.1,0.1)
            
            cpts4 = []
            cpts4.append([])
            cpts4.append([])
            cpts4.append([])
            cpts4[0].append(pt4)
            cpts4[0].append(pt3m)
            cpts4[1].append(pt4.mid(pt3))
            cpts4[1].append(pt3d)
            cpts4[2].append(pt3)
            cpts4[2].append(pt2m)
            ISurface(cpts4, 2, 1).clr(dratio/3+0.1,0.1,0.1)
            
            ISurface(pt1m, pt2, pt1).clr(dratio/3+0.1,0.1,0.1)
    surf.del()