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

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.

The problem here is for some reason, pygmsh does not support creating surface with more than 3 or 4 edges in a curve loop. However, the add_polygon function automatically creates the required surface and we can use that to build our 3D mesh.

The steps to create a 3D mesh from 2 2D cross-section are as follows:

1. Add the 2D cross sections as polygons to the mesh.
2. Reuse the points generated from the polygons for creating additional side surface lines.
3. Side surface lines are added as per suitable algorithm, here I use the minimum distance in the 2D projection between points to connect them with lines.
4. Surface lines along with cross section lines make up 2D surfaces, add each of the surfaces to the mesh.
5. Add a surface loop using the cross section and the side surfaces.
6. Add a volume from the surface loop.

The code snippet for the above steps is give below:

def calculate_distance_between_points(point1, point2):
    # calculate the projection distance between two points
    distance = ((point1.x[0] - point2.x[0]) ** 2 + (point1.x[1] - point2.x[1]) ** 2) ** 0.5
    return distance

def loft(polygon1, polygon2):
    with pygmsh.geo.Geometry() as geom:
        # set the larger polygon as polygon 1
        if len(polygon1.exterior.coords) < len(polygon2.exterior.coords):
            polygon1, polygon2 = polygon2, polygon1
        poly1 = geom.add_polygon(polygon1.exterior.coords[:-1])
        poly2 = geom.add_polygon(polygon2.exterior.coords[:-1])
    
        # for all the points in the first polygon, create a line to the corresponding point in the second polygon with the smallest distance
        lines = []
        lines_idx = []
        for idx_1, point1 in enumerate(poly1.points):
            min_distance = float("inf")
            for idx_2, point2 in enumerate(poly2.points):
                distance = calculate_distance_between_points(point1, point2)
                if distance < min_distance:
                    min_distance = distance
                    closest_point = poly2.points[idx_2]
                    idx_2_closest = idx_2
                lines.append([point1, closest_point])
                lines_idx.append([idx_1, idx_2_closest])

        side_surface_lines = [geom.add_line(p0 = poly1.points[line_idx[0]], p1 = poly2.points[line_idx[1]]) for line_idx in lines_idx]
    
        side_surfaces = []
        buffer_index = 0
        for idx, line in enumerate(poly1.lines):
            next_line_idx = (idx+1)%len(poly1.lines)
            # if the surface is rectangle create surface with 4 edges
            if side_surface_lines[next_line_idx].points[1] ==  side_surface_lines[idx].points[1]:
                curve_loop = geom.add_curve_loop([poly1.lines[idx], side_surface_lines[next_line_idx], -side_surface_lines[idx]])
                buffer_index -= 1
            else:
                curve_loop = geom.add_curve_loop([poly1.lines[idx], side_surface_lines[next_line_idx], -poly2.lines[idx+buffer_index], -side_surface_lines[idx]])
            side_surfaces.append(geom.add_surface(curve_loop))

        surface_loop = geom.add_surface_loop([poly1.surface] + side_surfaces + [poly2.surface])
        volume = geom.add_volume(surface_loop)

        mesh = geom.generate_mesh()
    return mesh

This creates the following 3D mesh: 3D mesh output