I’m new with gmsh and working with it to mesh a rocket nozzle. My teacher which is helping me with my thesis passed me this code, and the following error message is showing:

Exception: Surface 1 cannot be meshed using the transfinite algorithm (non-matching number of nodes on opposite sides 1524 != 1473 or 99 != 99)

Error : Surface 1 cannot be meshed using the transfinite algorithm (non-matching number of nodes on opposite sides 1524 != 1473 or 99 != 99)

I don’t have any idea what is going on! I just beginning to work with gmsh, and with python things get much more abstract. If anyone knows how to solve it, I’ll be eternaly greatful!

import gmsh
import sys
import numpy as np
import matplotlib.pyplot as plt

def neg(lst):
    return [ -i for i in lst ]

def getS(n,r):
    if r!=1:
        s=(r-1)/(r**n-1)
    elif r==1:
        s=1/n
    return s

def getPoints(start,stop,n,r):
    s=getS(n,r)
    x=np.zeros(n+1)
    x[0]=start
    x[1:]=start+(stop-start)*s*np.cumsum(np.array([r**i for i in range(n)]))
    return x
    
def chamber(Rc,step_count,Lc,progression):
    x = getPoints(0,Lc,step_count,progression)#np.linspace(0,Lc,step_count)
    y = Rc*np.ones(len(x))
    return [x,y]

def chamber_exit_circ(Rc,Rt,Lc,theta_t_i,step_count):
    chamber_e_l = Lc
    angles = np.linspace(90, 90 - theta_t_i, step_count+1)
    Xc = chamber_e_l
    Yc = Rc -1.5*Rt
    x = 1.5*Rt*np.cos(angles*np.pi/180) + Xc   
    y = 1.5*Rt*np.sin(angles*np.pi/180) + Yc 
    return [x,y]

def throat_inlet(Rt,theta_t_i,step_count,progression):
    throat_i_l = 1.5*Rt*np.sin(theta_t_i*np.pi/180)
    angles = getPoints(270 - theta_t_i, 270, step_count+1,progression)
    Xc = throat_i_l
    Yc = 2.5*Rt;
    x = 1.5*Rt*np.cos(angles*np.pi/180) + Xc
    y = 1.5*Rt*np.sin(angles*np.pi/180) + Yc
    return [x,y]

def throat_exit(Rt,theta_N,step_count):
    throat_e_l = 0.382*Rt*np.sin(theta_N*np.pi/180)
    angles = np.linspace(270 , 270 + theta_N ,step_count+1)
    Xc = throat_e_l
    Yc = 1.382*Rt
    x = 0.382*Rt*np.cos(angles*np.pi/180) + Xc 
    y = 0.382*Rt*np.sin(angles*np.pi/180) + Yc 
    return[x,y]

def bell_curve(Rt,Re,theta_N,theta_e,Ln_ratio,step_count,progression):
    # bezler quadratic curve equation %
    expansion_ratio = (Re/Rt)**2
    Nx = 0.382*Rt*np.cos((theta_N - 90)*np.pi/180)
    Ny = 1.382*Rt + 0.382*Rt*np.sin((theta_N - 90)*np.pi/180)
    Ex = Ln_ratio*( (np.sqrt(expansion_ratio) -1)*Rt + 1.5*Rt*(1/np.cos(15*np.pi/180) -1))/np.tan(15*np.pi/180) # Ln
    Ey = Re

    m1 = np.tan(theta_N*np.pi/180)
    m2 = np.tan(theta_e*np.pi/180)

    C1 = Ny - m1*Nx
    C2 = Ey - m2*Ex

    Qx = (C2 - C1)/(m1-m2)
    Qy = (m1*C2 -m2*C1)/(m1-m2)

    x = np.zeros(step_count+1)
    y = np.zeros(step_count+1)
    i=0

    for t in getPoints(0,1,step_count,progression):#np.linspace(0,1,step_count):
        x[i] = (1-t)**2*Nx + 2*(1-t)*t*Qx + t**2*Ex 
        y[i] = (1-t)**2*Ny + 2*(1-t)*t*Qy + t**2*Ey 
        i=i+1
        
    return [x,y]
 
# Rocket Dimensions

Rc = 71   # [mm]
Rt = 35.5 # [mm]
Re = 275  # [mm]

Lc = 100 # [mm]
Ln_ratio = 0.8; # Ratio of conical nozzle 

theta_t_i = 30 # [deg]
theta_N = 30   # [deg]
theta_e = 7    # [deg]

step_count_c_e = 40 #20;
step_count_t_i = 60
step_count_t_e = 60 #40 
step_count_r = 100 #80
progression_i=1/1.1
progression_bell=1.01
progression_r=1/1.05


# throat exit section
[x_t_e,y_t_e]=throat_exit(Rt,theta_N,step_count_t_e)

# throat inlet section
deltaN_t_i=x_t_e[1]-x_t_e[0]
throat_i_l = 1.5*Rt*np.sin(theta_t_i*np.pi/180)
delta0_t_i=deltaN_t_i*progression_i-throat_i_l *(progression_i-1)
step_count_t_i=int(np.ceil(np.log(deltaN_t_i/delta0_t_i)/np.log(progression_i)))
[x_t_i,y_t_i]=throat_inlet(Rt,theta_t_i,step_count_t_i,progression_i)

# chamber exit section
[x_c_e,y_c_e]=chamber_exit_circ(Rc,Rt,Lc,theta_t_i,step_count_c_e)
delta0_c_e=x_c_e[1]-x_c_e[0]

# chamber exit to throat inlet line
deltaN_ce_ti=x_t_i[1]-x_t_i[0]
L_c_e_ti=np.abs(np.tan((90 + theta_t_i)*np.pi/180)*(y_c_e[-1]-y_t_i[0]))
delta0_ce_ti=x_c_e[-1]-x_c_e[-2]
progression_ce_ti=(L_c_e_ti-delta0_ce_ti)/(L_c_e_ti-deltaN_ce_ti)
step_count_ce_ti=int(np.ceil(np.log(deltaN_ce_ti/delta0_ce_ti)/np.log(progression_ce_ti)))
y_ce_ti = getPoints(y_c_e[-1],y_t_i[0],step_count_ce_ti,progression_ce_ti)
x_ce_ti = np.tan((90 + theta_t_i)*np.pi/180)*y_ce_ti

# chamber section 
step_count_c=int(np.ceil(np.log(Lc*(1/progression_i-1)/delta0_c_e+1)/np.log(1/progression_i)))-1
[x_c,y_c]=chamber(Rc,step_count_c,Lc,progression_i)
delta0_c=Lc*getS(step_count_c,progression_i)

## chamber exit section (translate)
x_c_e = x_c_e + x_c[-1] - x_c_e[0]

#translate chamber exit to throat inlet, throat inlet and throat exit section sections
x_ce_ti = x_ce_ti + x_c_e[-1] - x_ce_ti[0]
x_t_i = x_t_i + x_ce_ti[-1] - x_t_i[0]
x_t_e = x_t_e + x_t_i[-1] - x_t_e[0]

# bell section
expansion_ratio = (Re/Rt)**2
delta0_t_e=x_t_e[-1]-x_t_e[-2]
Ln=Ln_ratio*( (np.sqrt(expansion_ratio) -1)*Rt + 1.5*Rt*(1/np.cos(15*np.pi/180) -1))/np.tan(15*np.pi/180)
step_count_nozzle=int(np.ceil(np.log(Ln*(progression_bell-1)/delta0_t_e+1)/np.log(progression_bell)))-1
[x_bell,y_bell]=bell_curve(Rt,Re,theta_N,theta_e,Ln_ratio,step_count_nozzle,progression_bell)
#translate bell section
x_bell = x_bell + x_t_e[-1] - x_bell[0]

gmsh.initialize()
gmsh.model.add("RaoNozzle")

factory=gmsh.model.occ

lc = 1
#pos=np.loadtxt('/Users/jsalazar/matlab/RAO-Bell-Nozzle/nozzleContour.csv',delimiter=',', skiprows=1)
x=np.concatenate((x_c,x_c_e[1:],x_ce_ti[1:],x_t_i[1:],x_t_e[1:],x_bell[1:]))
y=np.concatenate((y_c,y_c_e[1:],y_ce_ti[1:],y_t_i[1:],y_t_e[1:],y_bell[1:]))

pos=np.concatenate((x,y))
pos=pos.reshape((len(x),2),order='F')
npts=len(pos)
erase=[]

for i in range(0,npts-1):
    if np.linalg.norm(pos[i]-pos[i+1])==0:
        erase.append(i)

pos=np.delete(pos,erase,0)
npts=len(pos)

np.savetxt('nozzleContour.csv', pos/1e3,delimiter=',', header='x [m] y [m]')

pAxis=[[]]*npts
pBell=[[]]*npts

lAxis=[[]]*(npts-1)
lBell=[[]]*(npts-1)
lRadial=[[]]*2


for i in range(0,npts):
    pBell[i]=factory.addPoint(pos[i,0],pos[i,1],0,lc)  
    pAxis[i]=factory.addPoint(pos[i,0],0,0,lc)  
        
for i in range(0,npts-1):
    lBell[i]=factory.addLine(pBell[i],pBell[i+1])
    lAxis[i]=factory.addLine(pAxis[i],pAxis[i+1])

j=0


for i in [0,npts-1]:
    print(i)    
    lRadial[j]=factory.addLine(pAxis[i],pBell[i])
    factory.synchronize()
    gmsh.model.mesh.setTransfiniteCurve(lRadial[j], step_count_r, meshType="Progression", coef=progression_r)
    j=j+1

    
lAxis.reverse()    
loop=lBell+neg([lRadial[1]])+neg(lAxis)+[lRadial[0]]

cLoop=[]
cLoop.append(factory.addCurveLoop(loop))
surface=[]
surface.append(factory.addPlaneSurface(cLoop))
factory.synchronize()
gmsh.model.mesh.setRecombine(2, surface[0])
factory.synchronize()
gmsh.model.mesh.setTransfiniteSurface(surface[0],cornerTags=[pAxis[0],pAxis[-1],pBell[0],pBell[-1]])
factory.synchronize()
v=factory.extrude([(2,surface[0])],0,0,1,numElements=[1],recombine=True)
factory.synchronize()
gmsh.model.addPhysicalGroup(2,[surface[0],v[0][1]], name="frontAndBack")
gmsh.model.addPhysicalGroup(2,[v[-1][1]], name="inlet")
gmsh.model.addPhysicalGroup(2,[x[1] for x in v[2:npts+1]], name="wall")
gmsh.model.addPhysicalGroup(2,[v[npts+1][1]], name="outlet")
gmsh.model.addPhysicalGroup(2,[x[1] for x in v[npts+2:2*npts+1]], name="axis")
gmsh.model.addPhysicalGroup(3,[v[1][1]], name="internal")
                  
factory.synchronize()
gmsh.model.mesh.generate(3)
gmsh.write("raoNozzle.msh2")

if '-nopopup' not in sys.argv:
    gmsh.fltk.run()

gmsh.finalize()

I tryed to change the loop definitions, read the gmsh tutorial on python, changed a few things, but nothing solved the problem. I redownloaded the original python file and simply copied and pasted in this question, since I changed so many things…