I'm trying to dive a domain (in this case is a surface in 4 segments), this is how the surface is created and how I create the segments and check if a triangle is inside a segment

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
from mpl_toolkits.mplot3d import Axes3D
from utilidades import d2r, mag, r2d, test_creation
from sun import sun3, sun_reflection, normal_distribution, ray_count
from scipy.spatial import Delaunay, ConvexHull
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from matplotlib.tri import Triangulation




# Definimos los parámetros a, b, c
a = 40  # Semieje en x
b = 35  # Semieje en y
c = 20 # Semieje en z
dx=11

theta = np.linspace(d2r(30), d2r(50), dx) #ultimo parametro indica la cantidad de puntos de malla
phi = np.linspace(d2r(50),d2r(130),dx)
theta, phi = np.meshgrid(theta, phi)

x_e = a * np.cos(phi) * np.sin(theta)
y_e = b  * np.sin(phi) * np.sin(theta)
z_e = c  * np.cos(theta) 

##Superficie cilindrica
R=20 #supongamos que cm
h=x_e.min() #supongamos que cm

thet=np.linspace(d2r(0), d2r(75), dx)
x=np.linspace(-h, h, dx)
thet, x= np.meshgrid(thet, x)
y = R * np.cos(thet)
z = R * np.sin(thet)


# Aplanar los datos para formar una lista de puntos
points_2D = np.column_stack((theta.ravel(), phi.ravel()))  # Coordenadas paramétricas
points_3D = np.column_stack((x_e.ravel(), y_e.ravel(), z_e.ravel()))  # Coordenadas en la elipsoide


# Triangulación de Delaunay en el espacio paramétrico
tri = Delaunay(points_2D)

centr_tri = np.array([np.mean(points_3D[simplex], axis=0) for simplex in tri.simplices])

# Dividir el dominio paramétrico en 4 secciones
theta_mid = (theta.min() + theta.max()) / 2
phi_mid = (phi.min() + phi.max()) / 2

def create_division(first_parameter, second_parameter, mesh_poitns_2D, num_sections =4):
    middle_point_1 = (first_parameter.min() + first_parameter.max()) / (num_sections / 2)
    middle_point_2 = (second_parameter.min() + second_parameter.max()) / (num_sections / 2)


    segments = [
        mesh_poitns_2D[(first_parameter.ravel() <= middle_point_1) & (second_parameter.ravel() <= middle_point_2)],
        mesh_poitns_2D[(first_parameter.ravel() <= middle_point_1) & (second_parameter.ravel() >= middle_point_2)],
        mesh_poitns_2D[(first_parameter.ravel() >= middle_point_1) & (second_parameter.ravel() <= middle_point_2)],
        mesh_poitns_2D[(first_parameter.ravel() >= middle_point_1) & (second_parameter.ravel() >= middle_point_2)]
    ]

    return segments

def tri_in_segment(tri, sections):
    triangles =[]

    for i in range(len(sections)):
        tri_petenence = tri.find_simplex(sections[i])
        print(tri_petenence)
        triangles.append(tri_petenence)

    return triangles

seg = create_division(theta, phi, points_2D)

triangles = tri_in_segment(tri, seg)

When I plot the triangles, I can see that despite some triangles being inside the segment they are treated as if they aren't inside the segment

For better viewing I will show the polar representation of the surface