[edan0314]

I'm working with it to mesh a rocket nozzle. My teacher who 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 why it is not working, I spent my last 2 weeks in this problem. I'm just beginning to work with gmsh, and with python things get much more abstract. If anyone knows how to solve it, I'll be eternally grateful.

The real GMSH code begins at line 190, before that, it is just the definition of the nozzle geometry, and it is all correct.

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()