def _compute_VesPoly(R=2.4, r=1., elong=0., Dshape=0.,
                    divlow=True, divup=True, nP=200):
    """ Utility to compute three 2D (R,Z) polygons

    One represents a vacuum vessel, one an outer bumper, one a baffle

    The vessel polygon is centered on (R,0.), with minor radius r
    It can have a vertical (>0) or horizontal(<0) elongation in [-1;1]
    It can be D-shaped (Dshape in [0.,1.], typically 0.2)
    It can be non-convex, with:
        * a lower divertor-like shape
        * a upper divertor-like shape
    The elongation also affects the outer bumper and baffle

    Parameters
    ----------
    R:          int / float
        Major radius used as a center of the vessel
    r :         int / float
        Minor radius of the vessel
    elong:      int / float
        Dimensionless elongation parameter in [-1;1]
    Dshape:     int / float
        Dimensionless parameter for the D-shape (in-out asymmetry) in [0;1]
    divlow:     bool
        Flag indicating whether to incude a lower divertor-like shape
    divup:      bool
        Flag indicating whether to incude an upper divertor-like shape
    nP :        int
        Parameter specifying approximately the number of points of the vessel

    Return
    ------
    poly:       np.ndarray
        Closed (2,nP) polygon of the vacuum vessel, optionnally with divertors
    pbump:      np.ndarray
        Closed (2,N) polygon defining the outer bumper
    pbaffle:    np.ndarray
        Closed (2,N) polygon defining the lower baffle
    """

    # Basics (center, theta, unit vectors)
    cent = np.r_[R,0.]
    theta = np.linspace(-np.pi,np.pi,nP)
    poly = np.array([np.cos(theta), np.sin(theta)])

    # Divertors
    pdivR = np.r_[-0.1,0.,0.1]
    pdivZ = np.r_[-0.1,0.,-0.1]
    if divlow:
        ind = (np.sin(theta)<-0.85).nonzero()[0]
        pinsert = np.array([pdivR, -1.+pdivZ])
        poly = np.concatenate((poly[:,:ind[0]], pinsert, poly[:,ind[-1]+1:]),
                              axis=1)

    if divup:
        theta = np.arctan2(poly[1,:], poly[0,:])
        ind = (np.sin(theta)>0.85).nonzero()[0]
        pinsert = np.array([pdivR[::-1], 1.-pdivZ])
        poly = np.concatenate((poly[:,:ind[0]], pinsert, poly[:,ind[-1]+1:]),
                              axis=1)

    # Modified radius (by elongation and Dshape)
    rbis = r*np.hypot(poly[0,:],poly[1,:])
    theta = np.arctan2(poly[1,:],poly[0,:])
    rbis = rbis*(1+elong*0.15*np.sin(2.*theta-np.pi/2.))
    if Dshape>0.:
        ind = np.cos(theta)<0.
        coef = 1 + Dshape*(np.sin(theta[ind])**2-1.)
        rbis[ind] = rbis[ind]*coef

    er = np.array([np.cos(theta), np.sin(theta)])
    poly = cent[:,np.newaxis] + rbis[np.newaxis,:]*er

    # Outer bumper
    Dbeta = 2.*np.pi/6.
    beta = np.linspace(-Dbeta/2.,Dbeta/2., 20)
    pbRin = 0.85*np.array([np.cos(beta), np.sin(beta)])
    pbRout = 0.95*np.array([np.cos(beta), np.sin(beta)])[:,::-1]
    pinsert = np.array([[0.95,1.05,1.05,0.95],
                        [0.05,0.05,-0.05,-0.05]])

    ind = (np.abs(pbRout[1,:])<0.05).nonzero()[0]
    pbump = (pbRin, pbRout[:,:ind[0]], pinsert,
             pbRout[:,ind[-1]+1:], pbRin[:,0:1])
    pbump = np.concatenate(pbump, axis=1)
    theta = np.arctan2(pbump[1,:],pbump[0,:])
    er = np.array([np.cos(theta), np.sin(theta)])
    rbis = r*(np.hypot(pbump[0,:],pbump[1,:])
              *(1.+elong*0.15*np.sin(2.*theta-np.pi/2.)))
    pbump = cent[:,np.newaxis] + rbis[np.newaxis,:]*er

    # Baffle
    offR, offZ = 0.1, -0.85
    wR, wZ = 0.2, 0.05
    pbaffle = np.array([offR + wR*np.r_[-1,1,1,-1,-1],
                        offZ + wZ*np.r_[1,1,-1,-1,1]])
    theta = np.arctan2(pbaffle[1,:],pbaffle[0,:])
    er = np.array([np.cos(theta), np.sin(theta)])
    rbis = r*(np.hypot(pbaffle[0,:],pbaffle[1,:])
              *(1.+elong*0.15*np.sin(2.*theta-np.pi/2.)))
    pbaffle = cent[:,np.newaxis] + rbis[np.newaxis,:]*er

    return poly, pbump, pbaffle