The following example shows rectangular panelization with gaps and subdivision techniques of orthogonal grid.
The sample input file used in the example is below.

surface3.3dm

First, the code below shows rectangular panelization with gaps.

add_library('igeo')

size(480,360,IG.GL)

IG.open("surface3.3dm")
surfs = IG.surfaces()
surf = surfs[0]

unum = 10
vnum = 20
uinc = 1./unum
vinc = 1./vnum
# gap ratio between rectangular panel
gap = 0.25

for i in range(unum) : 
    for j in range(vnum) :
        pt1 = surf.pt(i*uinc, j*vinc)
        pt2 = surf.pt((i+1-gap)*uinc, j*vinc)
        pt3 = surf.pt((i+1-gap)*uinc, (j+1-gap)*vinc)
        pt4 = surf.pt(i*uinc, (j+1-gap)*vinc)
        ISurface(pt1,pt2,pt3,pt4)

surf.del()

Next, the code below shows the use of IBox providing 8 points as corners of the box.

add_library('igeo')

size(480,360,IG.GL)

IG.open("surface3.3dm")
surfs = IG.surfaces()
surf = surfs[0]

unum = 20
vnum = 40
uinc = 1./unum
vinc = 1./vnum
gap = 0.25

for i in range(unum) :
    for j in range(vnum) :
        
        pt1 = surf.pt(i*uinc, j*vinc)
        pt2 = surf.pt((i+1-gap)*uinc, j*vinc)
        pt3 = surf.pt((i+1-gap)*uinc, (j+1-gap)*vinc)
        pt4 = surf.pt(i*uinc, (j+1-gap)*vinc)
        
        # bottom points
        bpt1 = pt1.dup()
        bpt2 = pt2.dup()
        bpt3 = pt3.dup()
        bpt4 = pt4.dup()
        
        # move bottom points to xy-plane
        bpt1.z(0)
        bpt2.z(0)
        bpt3.z(0)
        bpt4.z(0)
        
        # create a box randomly 
        if IRand.pct(80) :
            IBox(pt1,pt2,pt3,pt4,bpt1,bpt2,bpt3,bpt4).clr(1-j*vinc)

surf.del()

Finally, extra for loops for the subdivision inside existing for loops with random number.

add_library('igeo')

size(480,360,IG.GL)

IG.open("surface3.3dm")
surfs = IG.surfaces()
surf = surfs[0]

unum = 20
vnum = 40
uinc = 1./unum
vinc = 1./vnum
gap = 0.15

for i in range(unum) : 
    for j in range(vnum) : 
        
        # random subdivision
        udiv=IRand.getInt(1,2)
        vdiv=IRand.getInt(1,4)
        uinc2 = 1./udiv
        vinc2 = 1./vdiv
        
        for k in range(udiv) : 
            for l in range(vdiv) : 
                
                pt1 = surf.pt((i + k*uinc2)*uinc, \
                              (j + l*vinc2)*vinc)
                pt2 = surf.pt((i+(k+1)*uinc2-gap )*uinc, \
                              (j + l*vinc2)*vinc)
                pt3 = surf.pt((i+(k+1)*uinc2-gap)*uinc, \
                              (j+(l+1)*vinc2-gap)*vinc)
                pt4 = surf.pt((i+k*uinc2)*uinc, \
                              (j+(l+1)*vinc2-gap)*vinc)
                
                # bottom points
                bpt1 = pt1.dup()
                bpt2 = pt2.dup()
                bpt3 = pt3.dup()
                bpt4 = pt4.dup()
                
                # move bottom points to xy-plane
                bpt1.z(0)
                bpt2.z(0)
                bpt3.z(0)
                bpt4.z(0)
                
                # lower the top randomly 
                if IRand.pct(20) :
                    pt1.sub(0,0,3)
                    pt2.sub(0,0,3)
                    pt3.sub(0,0,3)
                    pt4.sub(0,0,3)
                
                # create a box randomly 
                if IRand.pct(80) :
                    IBox(pt1,pt2,pt3,pt4,bpt1,bpt2,bpt3,bpt4).clr(IRand.gray())

surf.del()