Generate a mesh from my polygon geometry to iterate FEM for geometry optimization?

Question

The 2D script below generates polygons in a box, which represent a cross-section of three cylindrical rings parallel to the z axis which will have voltages applied. I'll solve the Laplace equation in cylindrical coordinates using the Finite Element Method as described in the answer to my SciComp SE question Simple, easy to install and use Python FEM solver (and example) for 2D cylindrical Laplace equation.

I'll convert the potential to an electric field and ray trace an electron beam through the lens, evaluate the results, and then iteratively adjust the lens geometry and voltages to optimize.

In this loop I'll need to convert my geometry to a proper mesh that scikit-fem can use. (The inside of the polygons is solid metal so I only need to make the mesh for the region in the box but outside the polygons).

I see that scikit-fem has built-in mesh constructors but from the page (and from the examples) I don't see how to take my polygons (with points and sides defined) and produce even a preliminary mesh. I'm completely new to this, a bit of hand-holding will be greatly appreciated!

For now I don't care what mesh generator I use, except that it needs to be really easy and foolproof to install (hopefully as simple pip install, I'm not a developer proper) and easily callable from within a python script.

Question: What is a simple automated way to generate a mesh from my polygon geometry so I can apply FEM (scikit-learn) in an iterative way to optimize my design?

The main script

V = 10

box = Rect_boundary('box', 0, 16, 0, 50, potential='GND', is_box=True)
lens_1 = Rect_boundary('lens 1', 13, 14.5, 8, 15, bevel=True, potential=+V)
lens_2 = Rect_boundary('lens 2', 13, 14.5, 18, 32, bevel=True, potential=-V)
lens_3 = Rect_boundary('lens 3', 13, 14.5, 35, 42, bevel=True, potential=+V)

things = box, lens_1, lens_2, lens_3

plot_things(things)
summarize_things(things)

A plot:

A printout:

thing:  box
is_box:  True
bevel:  False
points:
[[ 0  0]
 [50  0]
 [50 16]
 [ 0 16]]
facets:
[[0 1]
 [1 2]
 [2 3]
 [3 0]]

thing:  lens 1
is_box:  False
bevel:  True
points:
[[ 8.5 13. ]
 [14.5 13. ]
 [15.  13.5]
 [15.  14. ]
 [14.5 14.5]
 [ 8.5 14.5]
 [ 8.  14. ]
 [ 8.  13.5]]
facets:
[[0 1]
 [1 2]
 [2 3]
 [3 4]
 [4 5]
 [5 6]
 [6 7]
 [7 0]]

thing:  lens 2
is_box:  False
bevel:  True
points:
[[18.5 13. ]
 [31.5 13. ]
 [32.  13.5]
 [32.  14. ]
 [31.5 14.5]
 [18.5 14.5]
 [18.  14. ]
 [18.  13.5]]
facets:
[[0 1]
 [1 2]
 [2 3]
 [3 4]
 [4 5]
 [5 6]
 [6 7]
 [7 0]]

thing:  lens 3
is_box:  False
bevel:  True
points:
[[35.5 13. ]
 [41.5 13. ]
 [42.  13.5]
 [42.  14. ]
 [41.5 14.5]
 [35.5 14.5]
 [35.  14. ]
 [35.  13.5]]
facets:
[[0 1]
 [1 2]
 [2 3]
 [3 4]
 [4 5]
 [5 6]
 [6 7]
 [7 0]]

The classes and functions:

import numpy as np
import matplotlib.pyplot as plt

class Rect_boundary():
    def __init__(self, name, r1, r2, z1, z2, bevel=False, potential='ground',
                 is_box=False):
        """r2 > r1 and z2 > z1 please!"""
        self.name = name
        self.bevel = bevel
        self.is_box = is_box
        self.update_shape(r1, r2, z1, z2)
        self.update_potential(potential)
    def update_shape(self, r1, r2, z1, z2):
        self.r1 = r1
        self.r2 = r2
        self.z1 = z1
        self.z2 = z2
        self.thickness = self.r2 - self.r1
        b = self.thickness / 3.
        self.bevel_width = b
        if self.bevel:
            self.points = np.array([(z1+b, r1), (z2-b, r1), (z2, r1+b),
                                    (z2, r2-b), (z2-b, r2), (z1+b, r2),
                                    (z1, r2-b), (z1, r1+b)])
            self.facets = np.array([(0, 1), (1, 2), (2, 3), (3, 4),
                                    (4, 5), (5, 6), (6, 7), (7, 0)])
        else:
            self.points = np.array([(z1, r1), (z2, r1), (z2, r2), (z1, r2)])
            self.facets = np.array([(0, 1), (1, 2), (2, 3), (3, 0)])
        self.zmin, self.rmin = self.points.min(axis=0)
        self.zmax, self.rmax = self.points.max(axis=0)
        self.r_center = 0.5 * (self.rmin + self.rmax)
        self.z_center = 0.5 * (self.zmin + self.zmax)
        self.thickness = self.rmax - self.rmin
        self.length = self.zmax - self.zmin
    def update_potential(self, potential):
        if isinstance(potential, str) and potential.lower().startswith(('gr', 'gn')):
            self.potential = 0.
            self.color = 'black'
        else:
            self.ground = False
            self.potential = float(potential)
            if self.potential > 0:
                self.color = 'red'
            else:
                self.color = 'blue'

def plot_things(things):
    fs = 14

    fig, ax = plt.subplots(1, 1)
    for thing in things:
        x, y = zip(*thing.points)
        ax.plot(x, y, '.k')
        for i, j in thing.facets:
            (x0, y0), (x1, y1) = thing.points[i], thing.points[j]
            ax.plot([x0, x1], [y0, y1], '-', color=thing.color)
        if not thing.is_box:
            ax.text(thing.z_center, thing.rmin-3, thing.name,
                    color=thing.color, fontsize=fs, ha='center')
    ax.set_aspect('equal')
    ax.set_xlabel('z', fontsize=14)
    ax.set_ylabel('r', fontsize=14)
    plt.show()

def summarize_things(things):

    for thing in things:
        print('thing: ', thing.name)
        print('is_box: ', thing.is_box)
        print('bevel: ', thing.bevel)
        print('points:')
        print(thing.points)
        print('facets:')
        print(thing.facets)
        print('')