Gauss legendre for 2D for arbitrary points in xi,eta

Question

So im trying to code a general way to find the numerical integration gauss points in 2D. But want to choose my own points. All i can find on google is symmetrical points 2x2 and 3x3 etc. I have done it in 1D up to 5 points but im not sure how to do this is 2D. my code:

# Gauss Legendre Integrasjon 1D
from sympy import *
xi,eta = symbols("xi,eta")
Function = xi**6

D = int(input("Hvor mange dimensjoner? "))

if D == 1:

Vekter = [2, 1 , 8/9, 5/9, (18+sqrt(30))/(36), (18-sqrt(30))/(36), 128/225, (322+13*sqrt(70))/900, (322-13*sqrt(70))/900]

xi_p = [0, sqrt(1/3), -sqrt(1/3), 0, sqrt(3/5), -sqrt(3/5), sqrt((3/7)-(2/7)*sqrt(6/5)), -sqrt((3/7)-(2/7)*sqrt(6/5)), sqrt((3/7)+(2/7)*sqrt(6/5)), -sqrt((3/7)+(2/7)*sqrt(6/5)), 0, 1/3* sqrt(5-2*sqrt(10/7)), -1/3* sqrt(5-2*sqrt(10/7)), 1/3* sqrt(5+2*sqrt(10/7)), -1/3* sqrt(5+2*sqrt(10/7))  ]


N = int(input("Hvor mange integrasjons punkt? "))

if N == 1:
    
    INT_NUM = Vekter[0]*Function.subs(xi,xi_p[0])

if N == 2:
    
    INT_NUM = Vekter[1]*Function.subs(xi,xi_p[1])+Vekter[1]*Function.subs(xi,xi_p[2])
    
if N == 3:
    
    INT_NUM = Vekter[2]*Function.subs(xi,xi_p[3]) + Vekter[3]*Function.subs(xi,xi_p[4]) + Vekter[3]*Function.subs(xi,xi_p[5])
    
if N == 4:
    
    INT_NUM = Vekter[4]*Function.subs(xi,xi_p[6]) + Vekter[4]*Function.subs(xi,xi_p[7]) + Vekter[5]*Function.subs(xi,xi_p[8]) + Vekter[5]*Function.subs(xi,xi_p[9])
    
if N == 5:
    
    INT_NUM = Vekter[6]*Function.subs(xi,xi_p[10]) + Vekter[7]*Function.subs(xi,xi_p[11]) + Vekter[7]*Function.subs(xi,xi_p[12]) + Vekter[8]*Function.subs(xi,xi_p[13]) + Vekter[8]*Function.subs(xi,xi_p[14])

if D == 2:
N_xi = int(input("Hvor mange integrasjons punkt i xi? "))
N_eta = int(input("Hvor mange integrasjons punkt i eta? "))

Vekter = [2, 1 , 8/9, 5/9, (18+sqrt(30))/(36), (18-sqrt(30))/(36), 128/225, (322+13*sqrt(70))/900, (322-13*sqrt(70))/900]

xi_p = [0, sqrt(1/3), -sqrt(1/3), 0, sqrt(3/5), -sqrt(3/5), sqrt((3/7)-(2/7)*sqrt(6/5)), -sqrt((3/7)-(2/7)*sqrt(6/5)), sqrt((3/7)+(2/7)*sqrt(6/5)), -sqrt((3/7)+(2/7)*sqrt(6/5)), 0, 1/3* sqrt(5-2*sqrt(10/7)), -1/3* sqrt(5-2*sqrt(10/7)), 1/3* sqrt(5+2*sqrt(10/7)), -1/3* sqrt(5+2*sqrt(10/7))  ]
eta_p = [0, sqrt(1/3), -sqrt(1/3), 0, sqrt(3/5), -sqrt(3/5), sqrt((3/7)-(2/7)*sqrt(6/5)), -sqrt((3/7)-(2/7)*sqrt(6/5)), sqrt((3/7)+(2/7)*sqrt(6/5)), -sqrt((3/7)+(2/7)*sqrt(6/5)), 0, 1/3* sqrt(5-2*sqrt(10/7)), -1/3* sqrt(5-2*sqrt(10/7)), 1/3* sqrt(5+2*sqrt(10/7)), -1/3* sqrt(5+2*sqrt(10/7))  ]

if N_xi == 2 and N_eta ==2:
    for i in range(0,2):
        for n in range(0,2):
            print("s")


INT_NUM.evalf()