pygmsh: how to create 3D mesh from 2D cross-sections

Question

I have two 2D cross sections defined as a list of points stored in shapely.Polygon and I want to create a 3D mesh by lofting the first cross-section (cs1) to the second cross-section (cs2). This is similar to extrude, but instead of a simply adding depth d to a 2D shape, I want to join two 2D shapes at a given depth d.

I could not find an inbuilt function in pygmsh to accomplish this. However, pygmsh offer a range of geometry primitives that I may be able to use to get what I want.

Note, that the polygons cs1 and cs2 could be arbitrary and may contain holes. I don't know how I would accomplish this, but right now I am trying to solve the simplest case where I am trying to create a 3D mesh from two square cross sections as shown below. The 3d shape is encapsulated by the two squares and the black bold lines.

I tried extruding my 2D shapes using pygmsh.geo.Geometry().extrude but that is not something I want to do here.

I am attempting to manually create the 3D volume using the following steps:

1. From the polygons, get curve_loop
2. Use curve_loop to create surface
3. create surface loop from the list of surface defining the shell and holes.
4. create volume from the surface loop of shell and holes

But I have stumbled at the 2nd step to create surface from curve loop:

# Create the 2D polygons
poly1 = geom.add_polygon(polygon1.exterior.coords)
poly2 = geom.add_polygon(polygon2.exterior.coords)

# add curve loop to surface 
surf1 = geom.add_surface(poly1.curve_loop)

But I am getting the error:

Exception: Wrong definition of surface 3: 5 borders instead of 3 or 4

I do not know how to solve it.