
\section*{Librerie utilizzate}


\begin{minted}[fontsize=\footnotesize]{python}
# Codes transaltion
import numpy as np
from scipy import integrate
from math import gamma
from scipy.optimize import nnls

# Numerical results
import time
import statistics
import matplotlib.pyplot as plt

\end{minted}

\section*{{\tt{LHDM}} e {\tt{DM}}}

\begin{minted}[fontsize=\footnotesize]{python}

def LHDM(C, d, options=None, verbose=0):

# Code author: Laura Rinaldi
# https://github.com/laura-rinaldi/Cheap_Splinegauss

# Lawson-Hanson algorithm accelerated by Deviation Maximization (DM).

#     Parameters:
#         C (numpy.ndarray): Least squares matrix.
#         d (numpy.ndarray): Right-hand side vector.
#         options (dict, optional): Optimization parameters:
#             - init (bool): Use ULS initialization if True.
#             - tol (float): Tolerance for projected residual.
#             - k (int): Max number of indices added to Passive set
#               each iteration.
#             - thres (float): Threshold for angle between columns (0 to 1).
#         verbose (int, optional): Verbosity level.

#     Returns:
#         x (numpy.ndarray): Nonnegative solution minimizing ||C*x - d||.
#         resnorm (float): Squared residual norm ||C*x - d||^2.
#         exitflag (int): Exit condition (1: success, 0: exceeded iteration).
#         outeriter (int): Number of outer iterations.
#
    if options is None:
        options = {}

    m, n = C.shape
    nZeros = np.zeros(n)
    wz = np.zeros(n)
    itmax = 2 * m

    # Initialize sets
    P = np.zeros(n, dtype=bool)
    Z = np.ones(n, dtype=bool)
    x = np.zeros(n)

    thres = options.get('thres', 0.2222)
    thres_w = options.get('thres_w', 0.8)
    k = options.get('k', max(1, m // 20))
    tol = options.get('tol',10*np.finfo(float).eps*np.linalg.norm(C, 1)*len(C))
    LHDMflag = k > 1

    if verbose:
        print(f"LHDM({k}){'with ULS initialization'if\
        options.get('init',False)else ''}")

    if LHDMflag:
        Cnorm = C / np.linalg.norm(C, axis=0)

    # Initialize residual and dual variables
    resid = d - C @ x
    w = C.T @ resid
    outeriter = 0
    totiter = 0

    while np.any(Z) and (np.any(w[Z] > tol) or np.any(x[P] < 0)) and totiter < itmax:
        outeriter += 1
        totiter += 1

        wz[P] = -np.inf
        wz[Z] = w[Z]

        if outeriter == 1 or not LHDMflag:
            t = np.argmax(wz)
        else:
            t = DM(Cnorm, wz, k, thres, thres_w)
            t = t[:min(len(t), m - np.sum(P))]

        addedP = np.shape(t)
        z = np.zeros_like(x)
        P[t] = True
        Z[t] = False
        z[P] = np.linalg.lstsq(C[:, P], d, rcond=None)[0]

        iter = 0
        removedP = 0

        while np.any(z[P] <= 0) and totiter < itmax:
            totiter += 1
            iter += 1
            Q = (z <= 0) & P
            alpha = np.min(x[Q] / (x[Q] - z[Q]))
            x = x + alpha * (z - x)
            t = np.where((np.abs(x) < tol) & P)[0]
            removedP += len(t)
            Z[(np.abs(x) < tol) & P] = True
            P = ~Z
            z[P] = np.linalg.lstsq(C[:, P], d, rcond=None)[0]

        x = z
        resid = d - C @ x
        w = C.T @ resid

    exitflag = 1 if outeriter < itmax else 0
    resnorm = np.dot(resid, resid)

    return x #, resnorm, exitflag, outeriter




def DM(Cnorm, wz, k, thres, thres_w):

# Code author: Laura Rinaldi
# https://github.com/laura-rinaldi/Cheap_Splinegauss

#     Deviation Maximization

#     :param Cnorm: 2D NumPy array
#     :param wz: 1D NumPy array
#     :param k: integer
#     :param thres: threshold value
#     :param thres_w: weight threshold multiplier
#     :return: list of indices

    wzI = np.sort(wz)[::-1]  # Sort in descending order
    I = np.argsort(wz)[::-1]  # Get sorted indices
    t = I[0]
    p = [t]

    thres_wloc = thres_w * wzI[0]
    C = np.where(wzI > thres_wloc)[0]  # Get indices of wzI above threshold

    n = C.shape[0]
    add = 1

    for i in range(1, n):
        c = C[i]
        max_dev = np.max(np.abs(Cnorm[:, I[c]].T @ Cnorm[:, p]))

        if max_dev < thres:
            p.insert(0, I[c])
            add += 1

        if add >= k:
            break

    return p
    
\end{minted}

\section*{{\tt{Chebpolys}}}
\begin{minted}[fontsize=\footnotesize]{python}
def chebpolys(deg,x1):

# Code author: Giovanni Traversin
# Release date: 05 Sept 2025

# computes the Chebyshev-Vandermonde matrix on the real line by recurrence

# INPUT:
# deg = maximum polynomial degree
# x = list or array of abscissas

# OUTPUT:
# T = Chebyshev-Vandermonde matrix at x,
# T(i,j+1)=T_j(x_i), j=0,...,deg

    x = np.array(x1)
    n = len(x)
    T= np.zeros((n,deg+1))
    t0 = np.ones(n)
    t1 = x[:]
    T[:,0] = t0
    T[:,1] = t1

    for i in range(2,deg+1):
      t2 = 2*x*t1 - t0
      T[:,i] = t2
      t0 = t1[:]
      t1 = t2[:]

    return T

\end{minted}

\section*{{\tt{Chebvand}}}

\begin{minted}[fontsize=\footnotesize]{python}
def chebvand(deg, x, rect = None):

# Code author: Giovanni Traversin
# Release date: 05 Sept 2025

# INPUT:
# deg = polynomial degree
# x = 2-column array or 2 arrays of the same length of point coordinates
# rect = 4-component vector such that the rectangle
# [rect(1),rect(2)] x [rect(3),rect(4)] contains X
# If you do not compile the input for rect, the program will use as rect the
# smaller rectangle containing all the points in X.

# OUTPUT:
# V = Chebyshev-Vandermonde matrix at x, graded lexic. order

  X = np.array(x)
  if rect is None:
    rect = [min(X[:,0]), max(X[:,0]), min(X[:,1]), max(X[:,1])]

  # couples with length less or equal to deg
  # graded lexicographical order
  j = np.arange(deg + 1)
  j1, j2 = np.meshgrid(j,j)
  jj = j1 + j2
  dim = (deg+1)*(deg+2)//2
  couples = np.zeros((dim, 2), dtype=int)

  for s in range(deg+1):
    good = np.argwhere(jj == s)
    a = s*(s+1)//2
    for i  in range(len(good)):
      couples[a+i] = [j1[good[i][1],good[i][0]], j2[good[i][1],good[i][0]]]

  # mapping the mesh in the square [-1,1]^2
  a = rect[0]; b = rect[1]; c = rect[2]; d = rect[3]
  # map = [(2*X[:,0]-b-a)/(b-a) , (2*X[:,1]-d-c)/(d-c)]

  # Chebyshev-Vandermonde matrix on the mesh
  T1 = chebpolys(deg, (2*X[:,0]-b-a)/(b-a))
  T2 = chebpolys(deg, (2*X[:,1]-d-c)/(d-c))

  return  T1[:, couples[:, 0]] * T2[:, couples[:, 1]] # = V

\end{minted}

\section*{{\tt{Comprexcub}}}


\begin{minted}[fontsize=\footnotesize]{python}
def comprexcub(deg, x, omega, pos):

 # Code author: Giovanni Traversin
 # Release date: 05 Sept 2025

  # compression of bivariate cubature formulas by Tchakaloff points
  # or approximate Fekete points
  # useful, for example, in node reduction of algebraic cubature formulas
  # see the web page: http://www.math.unipd.it/~marcov/Tchakaloff.html

  # by Federico Piazzon, Alvise Sommariva and Marco Vianello
  # , May 2016


  # INPUT:
  # deg: polynomial exactness degree
  # X: 2-column array of point coordinates
  # omega: 1-column array of weights
  # pos: NNLS for pos=1, QR with column pivoting for pos=0

  # OUTPUT:
  # pts: 2-column array of extracted points
  # w: 1-column array of corresponding weights (positive for pos=1)
  # momerr: moment reconstruction error

  X = np.array(x)
  if len(X[0])>2:
    X = X.T
  rect = [min(X[:,0]), max(X[:,0]), min(X[:,1]), max(X[:,1])]
  V=chebvand(deg,X,rect)
  Q, R = np.linalg.qr(V)
  Q = np.real(Q)
  orthmom = Q.T @ omega
  match pos:
    case 1:
      weights = LHDM(Q.T,orthmom)
    case _:
      weights = nnls(Q.T,orthmom)[0]
  ind = np.where(abs(weights)>0)
  w = weights[ind]
  return X[ind, :][0], w, np.linalg.norm( Q[ind,:][0].T @ w - orthmom)
\end{minted}

\section*{{\tt{Gauss}}}

\begin{minted}[fontsize=\footnotesize]{python}
def gauss(N, albet):

# Code author: Giovanni Traversin
# Release date: 05 Sept 2025

#    Gauss quadrature rule.
#
#    Given a weight function w encoded by the nx2 array albet of the
#    first n recurrence coefficients for the associated orthogonal
#    polynomials, the first column of ab containing the n alpha-
#    coefficients and the second column the n beta-coefficients,
#    the call xw=GAUSS(n,ab) generates the nodes and weights xw of
#    the n-point Gauss quadrature rule for the weight function w.
#    The nodes, in increasing order, are stored in the first
#    column, the n corresponding weights in the second column, of
#    the nx2 array xw.

  ab = np.array(albet)
  if len(ab[:,0]) < N:
      raise ValueError('Input array ab is too short')
  J = np.zeros((N, N))

  for i in range(N):
      J[i][i] = ab[i][0]
  for j in range(N-1):
      J[j+1][j] = np.sqrt(ab[j+1][1])
      J[j][j+1] = J[j+1][j]

  D, V = np.linalg.eigh(J)
  I = np.argsort(D)
  D = np.sort(D)
  V = V[:,I]
  
  return np.transpose((D , ab[0][1]*V[0,:]**2))
  
\end{minted}


\section*{{\tt{R\_Jacobi}}}

\begin{minted}[fontsize=\footnotesize]{python}
def r_jacobi(N,a=0,b=[]):

# Code author: Giovanni Traversin
# Release date: 05 Sept 2025
  
#  Recurrence coefficients for monic Jacobi polynomials.
#
#    ab=R_JACOBI(n,a,b) generates the first n recurrence
#    coefficients for monic Jacobi polynomials with parameters
#    a and b. These are orthogonal on [-1,1] relative to the
#    weight function w(t)=(1-t)^a(1+t)^b. The n alpha-coefficients
#    are stored in the first column, the n beta-coefficients in
#    the second column, of the nx2 array ab. The call ab=
#    R_JACOBI(n,a) is the same as ab=R_JACOBI(n,a,a) and
#    ab=R_JACOBI(n) the same as ab=R_JACOBI(n,0,0).
#
#    Supplied by Dirk Laurie, 6-22-1998; edited by Walter
#    Gautschi, 4-4-2002.

  if b == []:
    b = a
  if N<=0 or a<=-1 or b<=-1:
    raise ValueError('parameter(s) out of range')
  nu = (b-a)/(a+b+2)
  mu = (2**(a+b+1))*gamma(a+1)*gamma(b+1)/gamma(a+b+2)

  # Important note: scipy.special.gamma() takes arbitrary np.arrays as input.
  # math.gamma() requires float or single-element numpy arrays, which a
  # crippling limitation in many use cases. –
  # The function scipy.special.gamma allows complex numbers, even though
  # math.gamma doesn't.

  if N == 1:
    return np.array([nu,mu])
  N0 = int(N)
  n = np.array(range(1,N0))
  nab = 2*n+a+b
  A = np.append((b**2-a**2)*np.ones(N0-1)/(nab*(nab+2)),nu)
  n = n[1:]
  nab = nab[1:]
  B1 = 4*(a+1)*(b+1)/((a+b+2)**2*(a+b+3))
  B = 4*(n+a)*(n+b)*n*(n+a+b)/((nab**2)*(nab+1)*(nab-1))
  ab = [mu] + [B1]
  for i in B:
    ab = ab + [i]
  ab = [A] + [ab]


  return   np.transpose(ab)
\end{minted}


\section*{{\tt{Chebyshev}}}

\begin{minted}[fontsize=\footnotesize]{python}

def chebyshev(N, Mom, Abm=False):

# Code author: Giovanni Traversin
# Release date: 05 Sept 2025

#  Modified Chebyshev algorithm.

#  Given a weight function w encoded by its first 2n modified
#  moments, stored in the (row) vector mom, relative to monic
#  polynomials defined by the (2n-1)x2 array abm of their
#  recurrence coefficients, [ab,normsq]=CHEBYSHEV(n,mom,abm)
#  generates the array ab of the first n recurrence coefficients
#  of the orthogonal polynomials for the weight function w, and
#  the vector normsq of their squared norms. The n alpha-
#  coefficients are stored in the first column, the n beta-
#  coefficients in the second column, of the nx2 array ab. The
#  call [ab,normsq]=CHEBYSHEV(n,mom) does the same, but using the
#  classical Chebyshev algorithm. If n is larger than the sizes
#  of mom and abm warrant, then n is reduced accordingly.

  mom = np.array(Mom)
  normsq = []
  if N <= 0:
     raise ('N out of range')
  if N > len(mom)/2:
    N = len(mom)/2

  if type(Abm) == bool:
    abm = np.zeros((2,int(2*N-1)))
  else:
    abm = np.array(Abm)
  if len(abm[0]) == 2 and len(abm)!=2:
      abm = abm.T

  if N > (len(abm[0])+1)/2:
    N = (len(abm[0])+1)/2

  if N == 1:
    normsq = [mom[0]]

  N = int(N)
  ab =  np.array([[abm[0][0]+mom[1]/mom[0], mom[0]]]+[[0,0]]*(N-1))
  sig = np.zeros((N+1,2*N))
  sig[1] = mom[0:int(2*N)]

  #NOTE: abm is a matrix with two raws instead of a matrix with two
  # columns, therefore the indexes are reversed.

  for n in range(2,N+1):
    for m in range(n-1,(2*N-n+1)):
      sig[n][m] = sig[n-1][m+1]-(ab[n-2][0]-abm[0][m])*sig[n-1][m]-\
      ab[n-2][1]*sig[n-2][m]+abm[1][m]*sig[n-1][m-1]
    ab[n-1,0] = abm[0][n-1] + sig[n][n]/sig[n][n-1] - \
    sig[n-1][n-1]/sig[n-1][n-2]
    ab[n-1,1] = sig[n][n-1]/sig[n-1][n-2]

  for i in range(N):
    normsq.append(sig[i+1][i])


  return ab, np.array(normsq)
\end{minted}


\section*{{\tt{Tridisolve}}}

\begin{minted}[fontsize=\footnotesize]{python}
def tridisolve(a,b,c,d):

  # Code author: Giovanni Traversin
  # Release date: 05 Sept 2025

  #   TRIDISOLVE  Solve tridiagonal system of equations.
  # From Cleve Moler's Matlab suite
  # http://www.mathworks.it/moler/ncmfilelist.html
  #     x = TRIDISOLVE(a,b,c,d) solves the system of linear equations
  #     b(1)*x(1) + c(1)*x(2) = d(1),
  #     a(j-1)*x(j-1) + b(j)*x(j) + c(j)*x(j+1) = d(j), j = 2:n-1,
  #     a(n-1)*x(n-1) + b(n)*x(n) = d(n).
  #
  #   The algorithm does not use pivoting, so the results might
  #   be inaccurate if abs(b) is much smaller than abs(a)+abs(c).
  #   More robust, but slower, alternatives with pivoting are:
  #     x = T\d where T = diag(a,-1) + diag(b,0) + diag(c,1)
  #     x = S\d where S = spdiags([[a; 0] b [0; c]],[-1 0 1],n,n)

  x = d
  n = len(x)-1
  for j in range(0,n):
    mu = a[j]/b[j]
    b[j+1] = b[j+1] - mu*c[j]
    x[j+1] = x[j+1] - mu*x[j]

  x[n] = x[n]/b[n]
  for i in range(n)[::-1]:
    x[i] = (x[i]-c[i]*x[i+1])/b[i]

  return x
\end{minted}


\section*{{\tt{Trigauss}}}

\begin{minted}[fontsize=\footnotesize]{python}
def trigauss(n, alpha, beta):

  # Code author: Giovanni Traversin
  # Release date: 05 Sept 2025
  
  # Original code by Gaspare Da Fies and Marco Vianello,
  # 8 Nov 2011

  # computes the n+1 angles and weights of a trigonometric gaussian
  # quadrature formula on [alpha,beta], 0<beta-alpha<=pi

  # uses the routines chebyshev.m, gauss.m from
  # W. Gautschi repository (not available online)
  # we suggest to put the following statements
  # ab = zeros(N,2); sig = zeros(N+1,2*N);
  # at the beginning of the body of chebyshev.m to speed-up execution

  # input:
  # n: trigonometric degree of exactness
  # [alpha,beta]: angular interval, 0<beta-alpha<=pi

  # output:
  # tw: (n+1) x 2 array of (angles,weights)

  # the formula integrates the canonical trigonometric basis with accuracy
  # from about 10^(-15) (small omega) to about 10^(-13) (omega-->pi)
  # up to n=300

  def integrand(t, n, omega):
    return np.cos(2*n*np.arccos(np.sin(t/2)/np.sin(omega/2)))

  omega=(beta-alpha)/2       # half-length of the angular interval

  #modified Chebyshev moments by recurrence
  z = np.zeros(n+1)
  z[0] = 2 * omega
  z[n] = integrate.quad(integrand,-omega, \
                         omega, args=(n, omega), limit=5000)[0]

  temp = np.arange(2, 2*n, 2)  # 2,4,...,2n-2
  dl = 1/4 - 1/(4*(temp-1))
  dc = 1/2 - 1/(np.sin(omega/2)**2) - 1/(2*(temp**2-1))
  du = 1/4 + 1/(4*(temp+1))
  d = 4 * np.cos(omega/2)/np.sin(omega/2)/(temp**2 - 1)
  d[n-2] = d[n-2] - du[n-2] * z[n]
  z[1:n] = tridisolve(dl[1:n-1], dc[0:n-1], du[0:n-2], d[0:n-1])
  mom = np.zeros(2*n+2)
  mom[0:2*n+1:2] = z  # fill even indices (0,2,4,...) with z

  # normalization of the moments (monic polynomials)
  k = np.arange(3, len(mom)+1)
  mom[2:] *= np.exp((2-k)*np.log(2))

  # recurrence coeffs of the monic Chebyshev polynomials
  abm = np.zeros((2*n+1, 2))
  abm[:,1] = 0.25
  abm[0,1] = np.pi
  abm[1,1] = 0.5

  # recurrence coeffs for the monic OPS w.r.t. the weight function
  # w(x)=2*sin(omega/2)/sqrt(1-sin^2(omega/2)*x^2) by the
  # modified Chebyshev algorithm
  ab, normsq = chebyshev(n+1,mom,abm)

  # Gaussian formula for the weight function above
  xw = gauss(n+1,ab)

  # angles and weights for the trigonometric gaussian formula
  return  np.column_stack((2*np.arcsin(np.sin(omega/2)*xw[:,0]) + \
   (beta+alpha)/2, xw[:,1]))
\end{minted}


\section*{{\tt{Gqellblend}}}

\begin{minted}[fontsize=\footnotesize]{python}
def gqellblend(n, A, B, C, alpha, beta):

  # Code author: Giovanni Traversin
  # Release date: 05 Sept 2025

  # Original code by Gaspare Da Fies, Alvise Sommariva and Marco Vianello

  # 2 June 2011

  # computes the nodes and weights of a product gaussian formula
  # exact on total-degree bivariate polynomials of degree <=n
  # on the planar region R obtained by linear blending (convex combination)
  # of two trigonometric arcs with parametric equations
  # P(theta)=A1*cos(theta)+B1*sin(theta)+C1
  # Q(theta)=A2*cos(theta)+B2*sin(theta)+C2
  # namely
  # R = {(x,y)=t*P(theta)+(1-t)*Q(theta), t in [0,1], theta in [alpha,beta],
  # 0<beta-alpha<=2*pi}

  # uses the routines:
  #
  # r_jacobi.m, gauss.m from
  # www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html
  #
  # trigauss.m
  # http://www.math.unipd.it/~marcov/mysoft/trigauss.m
  # this will be soon substituted by an optimized version input

  # input:
  # n: algebraic degree of exactness
  # A,B,C: 2x2 matrices of the parametric arc coefficients:
  # [alpha,beta]: angular interval, 0<beta-alpha<=2*pi

  # output:
  # xyw: 3 columns array of (xnodes,ynodes,weights)

  S1 = abs((A[0][0]-A[1][0])*(B[0][1]-B[1][1])+(A[0][1]-A[1][1])*\
   (B[1][0]-B[0][0]))>10*np.finfo(float).eps
  S2 = abs((C[0][0]-C[1][0])*(B[0][1]-B[1][1])+(C[0][1]-C[1][1])*\
   (B[1][0]-B[0][0]))>10*np.finfo(float).eps
  S3 = abs((A[0][0]-A[1][0])*(C[0][1]-C[1][1])+(A[0][1]-A[1][1])*\
   (C[1][0]-C[0][0]))>10*np.finfo(float).eps

  if S1 or S2 or S3:
    h=1
  else:
    h=0

  S4 = abs(A[0][1]*A[1][0]-A[0][0]*A[1][1]-B[0][1]*B[1][0]+\
           B[0][0]*B[1][1])>10*np.finfo(float).eps
  S5 = abs(A[0][1]*B[1][0]-A[0][0]*B[1][1]+B[0][1]*A[1][0]-\
           B[0][0]*A[1][1])>10*np.finfo(float).eps
  S6 = abs(B[1][0]*(C[0][1]-C[1][1])-B[1][1]*(C[0][0]-C[1][0]))\
  >10*np.finfo(float).eps
  S7 = abs(A[1][0]*(C[0][1]-C[1][1])-A[1][1]*(C[0][0]-C[1][0]))\
  >10*np.finfo(float).eps
  S8 = abs((C[0][0]-C[1][0])*(B[0][1]-B[1][1])+(C[0][1]-C[1][1])*\
   (B[1][0]-B[0][0]))>10*np.finfo(float).eps
  S9 = abs((A[0][0]-A[1][0])*(C[0][1]-C[1][1])+(A[0][1]-A[1][1])*\
   (C[1][0]-C[0][0]))>10*np.finfo(float).eps

  if S4 or S5:
    k = 2
  elif S6 or S7 or S8 or S9:
    k = 1
  else:
    k = 0

  # trigonometric gaussian formula on the arc
  tw=trigauss(n+k,alpha,beta)

  # algebraic gaussian formula on [0,1]
  ab=r_jacobi(int(np.ceil((n+h+1)/2)),0,0)

  xw=gauss(int(np.ceil((n+h+1)/2)),ab)
  xw[:,0]=xw[:,0]/2+1/2
  xw[:,1]=xw[:,1]/2

  #creating the grid
  t, theta = np.meshgrid(xw[:,0], tw[:,0])
  w1, w2 = np.meshgrid(xw[:,1], tw[:,1])

  # nodal cartesian coordinates and weights
  theta = theta.T.flatten()
  t = t.T.flatten()
  w1 = w1.T.flatten()
  w2 = w2.T.flatten()
  s = np.sin(theta)
  c = np.cos(theta)
  p1 = A[0][0]*c+B[0][0]*s+C[0][0]
  p2 = A[0][1]*c+B[0][1]*s+C[0][1]
  q1 = A[1][0]*c+B[1][0]*s+C[1][0]
  q2 = A[1][1]*c+B[1][1]*s+C[1][1]
  dp1 = -A[0][0]*s+B[0][0]*c
  dp2 = -A[0][1]*s+B[0][1]*c
  dq1 = -A[1][0]*s+B[1][0]*c
  dq2 = -A[1][1]*s+B[1][1]*c

  # xyw[:,0] = p1*t+q1*(1-t), p2*t+q2*(1-t)
  # xyw[:,1] = p2*t+q2*(1-t)
  # xyw[:,2] = np.abs((p1-q1)*(dp2*t+dq2*(1-t))-\
  # (p2-q2)*(dp1*t+dq1*(1-t)))*w1*w2)

  return  np.transpose((p1*t+q1*(1-t), p2*t+q2*(1-t), np.abs((p1-q1)*\
   (dp2*t+dq2*(1-t))-(p2-q2)*(dp1*t+dq1*(1-t)))*w1*w2))
\end{minted}


\section*{{\tt{Points2distances}}}



\begin{minted}[fontsize=\footnotesize]{python}

def points2distances(p):

  # Code author: Giovanni Traversin
  # Release date: 05 Sept 2025

  # Create euclidean distance matrix from point matrix.
  
  points = np.array(p)

  # All inner products between points.
  distances = points @ np.transpose(points)

  # Vector of squares of norms of points.
  lsq = np.diag(distances)

  a = np.array([[i]*len(points) for i in lsq])

  # distances = np.sqrt(a+np.transpose(a)-2*distances)
  return np.sqrt(a+np.transpose(a)-2*distances)
\end{minted}


\section*{{\tt{Cubature\_rules\_1D}}}



\begin{minted}[fontsize=\footnotesize]{python}

def cubature_rules_1D(n,cubature_type):

  # Code author: Giovanni Traversin
  # Release date: 05 Sept 2025

  # SEE WALDVOGEL PAPER. ADDED NODES

  # Weights of the Fejer2, Clenshaw-Curtis and Fejer1 quadrature by DFTs
  # n>1. Nodes: x_k = cos(k*pi/n)

  N = np.array(range(1,n,2))
  l = len(N)
  K = np.array(range(0,n-l))

  match cubature_type:
    case 1:       # FEJER 1
      v0 = np.append([2*np.exp(1j*np.pi*K/n)/(1-4*K**2)],[0]*(l+1))
      v1 = v0[:-1]+np.conjugate(v0[:0:-1])
      weights = np.real(np.fft.ifft(v1))  # Imaginary residual
      k = np.arange(.5 , n+.5)
      nodes = np.cos(k*np.pi/n)

    case 2:       # FEJER 2
      v0 = np.append(2/N/(N-2), [1/N[-1]]+ [0]*(n-l))
      v2 = -v0[:-1]-v0[:0:-1]
      weights = np.append(np.fft.ifft(v2),0)
      k = np.array(range(0,n+1))
      nodes = np.cos(k*np.pi/n)

    case 3:       # CLENSHAW CURTIS
      g0 = -np.ones(n)
      g0[l] = g0[l]+n
      g0[n-l] = g0[n-l]+n
      g=g0/(n**2-1+n%2)
      v0 = np.append(2/N/(N-2), [1/N[-1]]+ [0]*(n-l))
      v2 = -v0[:-1]-v0[:0:-1]
      wcc = np.real(np.fft.ifft(v2+g))  # Imaginary residual
      weights = np.append(wcc, wcc[0])
      k = np.array(range(0,n+1))
      nodes = np.cos(k*np.pi/n)

    case 4:       # GAUSS LEGENDRE
      beta = .5/np.sqrt(1-1/((2*np.array(range(1,n+1)))**2))
      T = np.diag(beta,1)+np.diag(beta,-1)
      nodes,V=np.linalg.eigh(T)
      index = np.argsort(nodes)
      nodes = np.sort(nodes)
      weights = 2*V[0]**2


  return nodes, weights
\end{minted}

\section*{{\tt{auto\_rotation}}}



\begin{minted}[fontsize=\footnotesize]{python}

def auto_rotation(p, V1, V2):

  # Code author: Giovanni Traversin
  # Release date: 05 Sept 2025

  # AUTOMATIC ROTATION OF A CONVEX POLYGON SO THAT "GAUSSIAN POINTS",
  # AS IN THE PAPER THEY ARE ALL CONTAINED IN THE CONVEX POLYGON.
  # SEE THE PAPER FOR DETAILS.

  polygon_bd = np.array(p)
  vertex_1 = np.array(V1)
  vertex_2 = np.array(V2)

  # FIND DIRECTION AND ROTATION ANGLE.
  if len(vertex_1) == 0:
    # COMPUTING ALL THE DISTANCES BETWEEN POINTS.
    # A LITTLE TIME CONSUMING AS PROCEDURE.
    distances = points2distances(polygon_bd)
    max_distances = distances.max(1)
    max_distance = distances.max()
    max_row_comp = np.argmax(max_distances)
    max_col_comp = []
    for i in distances:
      max_col_comp.append(np.argmax(i))
    vertex_1 = np.array(polygon_bd[max_col_comp[max_row_comp]])
    vertex_2 = np.array(polygon_bd[max_row_comp])
    direction_axis = (vertex_2-vertex_1)/max_distance
  else:
    direction_axis = (vertex_2-vertex_1)/np.linalg.norm(vertex_2-vertex_1)

  rot_angle_x = np.arccos(direction_axis[0])
  rot_angle_y = np.arccos(direction_axis[1])

  if rot_angle_y <= np.pi/2:
    if rot_angle_x <= np.pi/2:
      rot_angle = -rot_angle_y
    else:
      rot_angle = rot_angle_y
  else:
    if rot_angle_x <= np.pi/2:
      rot_angle = np.pi-rot_angle_y
    else:
      rot_angle = rot_angle_y

  # CLOCKWISE ROTATION.
  rot_matrix = np.array([[np.cos(rot_angle), np.sin(rot_angle)], \
   [-np.sin(rot_angle), np.cos(rot_angle)]])
  number_sides = len(polygon_bd[:,0])-1
  polygon_bd_rot = (rot_matrix @ polygon_bd.T).T
  axis_abscissa = rot_matrix @ vertex_1.T

  return polygon_bd_rot, rot_matrix, rot_angle, axis_abscissa, \
  vertex_1, vertex_2
\end{minted}

\section*{{\tt{Polygauss}}}

\begin{minted}[fontsize=\footnotesize]{python}


def polygauss(N, p, P1=[], Q1=[], rotation = 1):

  # Code author: Giovanni Traversin
  # Release date: 05 Sept 2025

  #--------------------------------------------------------------------------
  # REFERENCE PAPER:
  # [1] A. SOMMARIVA and M. VIANELLO
  # "Gauss-like and triangulation-free cubature over polygons".
  #
  # INPUT:
  #
  # N     : DEGREE OF THE 1 DIMENSIONAL GAUSS-LEGENDRE RULE.
  #
  # polygon_sides: IF THE POLYGON HAS "L" SIDES, "boundary.pts" IS A
  #         VARIABLE CONTAINING ITS VERTICES, ORDERED COUNTERCLOCKWISE.
  #         AS LAST ROW MUST HAVE THE COMPONENTS OF THE FIRST VERTEX.
  #         IN OTHER WORDS, THE FIRST ROW AND LAST ROW ARE EQUAL.
  #         "polygon_sides" IS A "L+1 x 2" MATRIX.
  #
  #            --------- NOT MANDATORY VARIABLES ---------
  #
  # rotation: 0: NO ROTATION.
  #           1: AUTOMATIC.
  #           2: PREFERRED DIRECTION ROTATION BY P, Q.
  #
  # P, Q: DIRECTION THAT FIXES THE ROTATION.
  #
  # OUTPUT:
  #
  # xyw     : THE GAUSS LIKE FORMULA PRODUCES THE NODES (xyw(:,1),xyw(:,2))
  #           AND THE WEIGHTS xyw(:,3) OF A CUBATURE RULE ON THE POLYGON.
  #
  #--------------------------------------------------------------------------
  # EXAMPLE 1 (NO ROTATION.)
  #---------------------------
  #
  # >> xyw=polygauss_2013(2,[0 0; 1 0; 1 1; 0 1; 0 0],0)
  #
  # xyw =
  #
  #     0.2113    0.2113    0.2500
  #     0.2113    0.7887    0.2500
  #     0.7887    0.2113    0.2500
  #     0.7887    0.7887    0.2500
  #
  # >>
  #
  #--------------------------------------------------------------------------
  # EXAMPLE 2 (AUTO ROTATION.)
  #-----------------------------
  #
  # >> xyw=polygauss_2013(2,[0 0; 1 0; 1 1; 0 1; 0 0])
  #
  # xyw =
  #
  #     0.0683    0.0444    0.0078
  #     0.3028    0.1972    0.0556
  #     0.5374    0.3499    0.0616
  #     0.6501    0.4626    0.0616
  #     0.8028    0.6972    0.0556
  #     0.9556    0.9317    0.0078
  #     0.9317    0.9556    0.0078
  #     0.6972    0.8028    0.0556
  #     0.4626    0.6501    0.0616
  #     0.3499    0.5374    0.0616
  #     0.1972    0.3028    0.0556
  #     0.0444    0.0683    0.0078
  #     0.1008    0.0119    0.0078
  #     0.4472    0.0528    0.0556
  #     0.7935    0.0938    0.0616
  #     0.9062    0.2065    0.0616
  #     0.9472    0.5528    0.0556
  #     0.9881    0.8992    0.0078
  #     0.8992    0.9881    0.0078
  #     0.5528    0.9472    0.0556
  #     0.2065    0.9062    0.0616
  #     0.0938    0.7935    0.0616
  #     0.0528    0.4472    0.0556
  #     0.0119    0.1008    0.0078
  #
  # >>

  #--------------------------------------------------------------------------
  # Copyright (C) 2007-2013 Marco Vianello and Alvise Sommariva
  #
  # This program is free software; you can redistribute it and/or modify
  # it under the terms of the GNU General Public License as published by
  # the Free Software Foundation; either version 2 of the License, or
  # (at your option) any later version.
  #
  # This program is distributed in the hope that it will be useful,
  # but WITHOUT ANY WARRANTY; without even the implied warranty of
  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  # GNU General Public License for more details.
  #
  # You should have received a copy of the GNU General Public License
  # along with this program; if not, write to the Free Software
  # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  #
  # Authors:
  # Marco Vianello    <marcov@euler.math.unipd.it>
  # Alvise Sommariva  <alvise@euler.math.unipd.it>
  # Date: April 30, 2013.
  #--------------------------------------------------------------------------

  polygon_sides = np.array(p)
  P = np.array(P1)
  Q = np.array(Q1)
  x_bd = polygon_sides[:,0]
  y_bd = polygon_sides[:,1]

  # "MINIMUM" RECTANGLE CONTAINING POLYGON.
  x_min = min(x_bd)
  x_max=max(x_bd)
  y_min = min(y_bd)
  y_max=max(y_bd);

  cubature_type=4

  # POLYGON ROTATION (IF NECESSARY).
  match rotation:
    case 0:
      #print('\n \t [ROTATION]: NO.')
      rot_matrix = np.eye(2)
      axis_abscissa = np.array([x_min, y_max])-np.array([x_min, y_min])
    case 1:
      #print('\n \t [ROTATION]: AUTOMATIC')
      polygon_sides, rot_matrix, rot_angle, axis_abscissa, P, Q = \
      auto_rotation(polygon_sides,[],[])
    case 2:
      #print('\n \t [ROTATION]: PREFERRED DIRECTION');
      nrm_vect = np.linalg.norm(Q-P)
      if nrm_vect > 0:
        direction_axis = (Q-P)/nrm_vect
        polygon_sides,rot_matrix,rot_angle,axis_abscissa,P,Q = \
        auto_rotation(polygon_sides,P,Q)
      else:
        #print('\n \t [WARNING]: THE DIRECTION VECTOR IS NULL. ')
        #print('USING AUTOMATIC ROTATION.')
        polygon_sides,rot_matrix,rot_angle,axis_abscissa,P,Q = \
        auto_rotation(polygon_sides,P,Q)


  # COMPUTE NODES AND WEIGHTS OF 1D GAUSS-LEGENDRE RULE.
  # TAKEN FROM TREFETHEN PAPER "Is ... Clenshaw-Curtis?".

  # DEGREE "N" (ORDER GAUSS PRIMITIVE)
  s_N, w_N = cubature_rules_1D((N-1),cubature_type)
  N_length = len(s_N)

  # DEGREE "M" (ORDER GAUSS INTEGRATION)
  M = N+1
  s_M, w_M = cubature_rules_1D((M-1),cubature_type)

  # L: NUMBER OF SIDES OF THE POLYGON.
  L = len(polygon_sides[:,1])-1

  #a=0.5;
  a = axis_abscissa[0]

  # COMPUTE 2D NODES (nodes_x,nodes_y) AND WEIGHTS "weights".
  nodes_x=[]
  nodes_y=[]
  weights=[]

  for i in range(L):
     x1=polygon_sides[i][0]
     x2=polygon_sides[i+1][0]
     y1=polygon_sides[i][1]
     y2=polygon_sides[i+1][1]
     if not (x1 == a and x2 == a):
      if y2-y1 != 0:
        if x2-x1 != 0:
          s_M_loc = s_M
          w_M_loc = w_M
        else:
          s_M_loc = s_N
          w_M_loc = w_N

        M_length = len(s_M_loc)
        half_pt_x = (x1+x2)/2
        half_pt_y = (y1+y2)/2
        half_length_x = (x2-x1)/2
        half_length_y = (y2-y1)/2

        #  GAUSSIAN POINTS ON THE SIDE.
        x_gauss_side = half_pt_x+half_length_x*s_M_loc; #SIZE: (M_loc,1)
        y_gauss_side = half_pt_y+half_length_y*s_M_loc; #SIZE: (M_loc,1)

        scaling_fact_plus = (x_gauss_side+a)/2  #SIZE: (M_loc,1)
        scaling_fact_minus = (x_gauss_side-a)/2  #SIZE: (M_loc,1)
        # SIZE:(M_loc,1)
        local_weights = (half_length_y*scaling_fact_minus)*w_M_loc
        # SIZE: (M_loc,N)
        term_1 = np.tile(np.transpose([scaling_fact_plus]),N_length)
        # SIZE: (M_loc,N)
        term_2 = np.tile(np.transpose([scaling_fact_minus]),N_length)
        rep_s_N = np.tile(s_N,(M_length,1))


        # x, y ARE STORED IN MATRICES. A COUPLE WITH THE SAME INDEX
        # IS A POINT, i.e. "P_i=(x(k),y(k))" FOR SOME "k".
        x= term_1+term_2*rep_s_N
        number_rows=len(x)
        number_cols=len(x[0])
        x = x.flatten('F')
        y = np.tile(y_gauss_side,N_length)
        # THE INVERSE OF A ROTATION MATRIX IS ITS TRANSPOSE.
        rot_gauss_pts = rot_matrix.T @ np.array([x,y])


        # GAUSS POINTS IN THE ORIGINAL SYSTEM.
        x_rot = rot_gauss_pts[0,:]
        y_rot = rot_gauss_pts[1,:]

        x_rot = np.reshape(x_rot,(number_cols,number_rows)).T
        y_rot = np.reshape(y_rot,(number_cols,number_rows)).T

        weights = np.append(weights, local_weights)
        for i in x_rot:
          nodes_x.append(i)
        for j in y_rot:
          nodes_y.append(j)

  #nodes_x = np.ravel(np.array(nodes_x),'F')
  #nodes_y = np.ravel(np.array(nodes_y),'F')
  #weights = np.ravel(np.outer(weights, w_N),'F')

  #return nodes_x, nodes_y, weights

  return np.transpose((np.ravel(np.array(nodes_x),'F'), \
                       np.ravel(np.array(nodes_y),'F'), \
                       np.ravel(np.outer(weights, w_N),'F')))
\end{minted}

\section*{{\tt{Circtrap}}}

\begin{minted}[fontsize=\footnotesize]{python}
def circtrap(n,a,b,c,d,cc,r):

  # Code author: Giovanni Traversin
  # Release date: 05 Sept 2025

  # by E. Artioli, A. Sommariva, M. Vianello
  # April 2018

  # n: polynomial exactness degree
  # a,b: circular arc extrema coords 1 x 2
  # c,d: base segment extrema coords 1 x 2, ac and bd are the sides
  # cc: circle center coords 1 x 2
  # r: circle radius

  # xyw: 3-column array xyw(:,1:2) nodes, xyw(:,3) weights 
  
  A = np.array(a)
  B = np.array(b)
  C = np.array(c)
  D = np.array(d)
  CC = np.array(cc)

  Z = (A[0]-CC[0])+1j*(A[1]-CC[1])
  W = (B[0]-CC[0])+1j*(B[1]-CC[1])
  az = np.angle(Z)
  aw = np.angle(W)

  if az<=aw:
    if aw-az<=np.pi:
      alpha = az
      beta = aw
      U = C
      V = D
    else:
      alpha = aw
      beta = az+2*np.pi
      U = D
      V = C

  if az>aw:
    if az-aw<=np.pi:
      alpha = aw
      beta = az
      U = D
      V = C
    else:
      alpha = az
      beta = aw+2*np.pi
      U = C
      V = D

  om = (beta-alpha)/2
  g = (beta+alpha)/2
  s=2*np.sin(om)
  A1 = [[r*np.cos(g), r*np.sin(g)], [0, 0]]
  B1 = [[-r*np.sin(g), r*np.cos(g)], [(V[0]-U[0])/s , (V[1]-U[1])/s ]]
  C1 = [[CC[0], CC[1]], [(V[0]+U[0])/2, (V[1]+U[1])/2]]

  return gqellblend(n,A1,B1,C1,-om,om)
\end{minted}

\section*{{\tt{Polygcirc}}}

\begin{minted}[fontsize=\footnotesize]{python}
def polygcirc(n,V,A1,B1,center,r,conv,pos=1):

  # Code author: Giovanni Traversin
  # Release date: 05 Sept 2025

  #--------------------------------------------------------------------------
  # OBJECT:

  # Computation of a basic and a compressed positive cubature formula
  # on a polygonal element with a circular side, namely the set
  # union (convex arc) or difference (concave arc) of a convex
  # polygonal element with a circular segment.
  #--------------------------------------------------------------------------
  # INPUT:

  # n: polynomial degree of exactness
  # A1,B1: extrema of the circular arc
  # center,r: circular arc center and radius
  # V: polygon vertices (2-column array of coords)
  #    note that to "v" are added by default the arc extrema, in
  #    couterclockwise order
  # conv: conv=1 for a convex arc, conv=0 for a concave arc
  #
  # WARNING: the figure vertices are a,v(1,:),...,v(end,:),b and MUST BE
  # in COUNTERCLOCKWISE order
  # (the arc ba is clockwise on the circle if concave and counterclockwise
  # if convex)
  #--------------------------------------------------------------------------
  # OUTPUT:

  # xyw: 3-column array of cubature nodes and positive weights
  # xywc: 3-column array of compressed cubature nodes and positive weights
  #--------------------------------------------------------------------------
  # AUTHORS:

  # Authors: E. Artioli, A. Sommariva and M. Vianello
  # Written: April 24, 2018
  # Revised: December 02, 2021
  #--------------------------------------------------------------------------

  # auxiliary function used to stack columns on xyw more easily

  cc = np.array(center)
  v = np.array(V)
  a = np.array(A1)
  b = np.array(B1)
  L = []
  P = []
  subs = []

  def stackxyw(xyw, nw):
    if len(xyw) == 0:
      return nw
    else:
      return np.vstack((xyw, nw))

  # Convex scenario
  if conv==1:
    P = np.vstack([a,v,b,a])
    xyw = np.vstack((polygauss(n,P), circtrap(n,b,a,b,a,cc,r)))
    L += [xyw[:,1].size]
    A = b[:]
    B = a[:]
    C = b[:]
    D = a[:]
    subs += [A,B,C,D]

  #Concave scenario
  else: #conv = 0
    k = 0; l = 0; t = 0; e=[]; z=[]; u=[]; eta=[]; xyw=[]
    aa=a-cc; bb=b-cc
    angleab = np.arccos((aa@bb)/(np.linalg.norm(aa)*np.linalg.norm(bb)))
    anglea = np.angle(aa[0]+1j*aa[1])
    angleb = np.angle(bb[0]+1j*bb[1])

    if angleb<=np.pi:
      clockba = anglea>=angleb and anglea<angleb+np.pi
    else:
      clockba = anglea>angleb-np.pi and anglea<=angleb
    if len(v[:,0])==2:
      v = np.array([v[0], (v[0]+v[1])/2, v[1]])

    vv = v-np.matlib.repmat(cc, len(v[:,0]), 1)

    for i in range(len(v[:,1])):
      angle1 = np.arccos((vv[i]@aa)/(np.linalg.norm(aa)*np.linalg.norm(vv[i])))
      angle2 = np.arccos((vv[i]@bb)/(np.linalg.norm(bb)*np.linalg.norm(vv[i])))

      if angle1<angle2 and angle2>angleab and clockba==0:
        e += [v[i]]
        k+=1
      if clockba==1 or (angle1<=angleab and angle2<=angleab):
        z += [v[i]]
        zz = z[l]-cc
        theta = np.angle(zz[0]+1j*zz[1])
        u += [cc+r*np.array([np.cos(theta), np.sin(theta)])]
        l+=1
      if angle2<angle1 and angle1>angleab and clockba==0:
        eta += [v[i]]
        t+=1

    e = np.array(e); z = np.array(z); u = np.array(u); eta  = np.array(eta)

    if k>1:
      P1 = np.vstack((a, e, a))
      xyw = polygauss(n, P1)
      L += [xyw[:,0].size]

    if k==0 and l>=1:
      if np.linalg.norm(u[0]-a)>10**(-14):
        nw = circtrap(n,a,u[0], z[0], z[0], cc, r)
        L += [nw[:,0].size]
        xyw = nw
        A = a; B = u[0,:]
        C = z[0]; D = z[0]
        subs += [A,B,C,D]

    if k>=1 and l>=1:
      nw = circtrap(n,a,u[0],e[k-1],z[0],cc,r)
      L += [nw[:,0].size]
      xyw = stackxyw(xyw, nw)
      A = a; B = u[0,:]
      C = e[k-1,:]; D = z[0,:]
      subs += [A,B,C,D]

    for j in range(1,l):
      nw=circtrap(n,u[j-1],u[j],z[j-1],z[j],cc,r)
      L += [nw[:,0].size]
      xyw = stackxyw(xyw, nw)
      A = u[j-1]; B = u[j]
      C = z[j-1]; D = z[j]
      subs += [A,B,C,D]

    if t == 0 and l>=1:
      if np.linalg.norm(u[l-1]-b)>10**(-14):
        nw = circtrap(n,u[l-1],b,z[l-1,],z[l-1],cc,r)
        L += [nw[:,0].size]
        xyw = stackxyw(xyw, nw)
        A = u[l-1]; B = b
        C = z[l-1]; D = z[l-1]
        subs += [A,B,C,D]

    if t>=1 and l>=1:
      nw = circtrap(n,u[l-1,:],b,z[l-1,:],eta[0,:],cc,r)
      L += [nw[:,0].size]
      xyw = stackxyw(xyw, nw)
      A = u[l-1]; B = b
      C = z[l-1]; D = eta[0]
      subs += [A,B,C,D]

    if t>1:
      P2=np.vstack((b,eta,b))
      nw=polygauss(n,P2)
      L += [nw[:,1].size]
      xyw = stackxyw(xyw, nw)

    if l==0:
      if k>=1 and t>=1:
        nw = circtrap(n,a,b,e[k-1],eta[0],cc,r)
        L += [nw[:,0].size]
        A = a; B = b
        C = e[k-1]; D = eta[0]
        subs += [A,B,C,D]

      if k==0 and t>=1:
        nw = circtrap(n,a,b,eta[0],eta[0],cc,r)
        L += [nw[:,0].size]
        A = a; B = b
        C = eta[0]; D = eta[0]
        subs += [A,B,C,D]

      if k>=1 and t==0:
        nw = circtrap(n,a,b,e[k-1],e[k-1],cc,r)
        L += [nw[:,0].size]
        A = a; B = b
        C = e[k-1,:]; D = e[k-1,:]
        subs += [A,B,C,D]
      xyw = stackxyw(xyw, nw)


  subs = np.array(subs)
  pts, w, momerr = comprexcub(n,[xyw[:,0], xyw[:,1]],xyw[:,2],pos)
  xywc = np.transpose([pts[:,0], pts[:,1], w])

  return  xyw, xywc, P, L, subs
\end{minted}

\section*{demo_polygcirc}
\begin{minted}[fontsize=\footnotesize]{python}
def demo_polygcirc(test_type=2,n=4,pos=2):

  # Code author: Giovanni Traversin
  # Release date: 05 Sept 2025

  # polynomial exactness degree
  # compression type: pos=0: lsqnonneg, pos=1: LHDM.
  # note: pos=0: may have negative weights
  # tests to determine median cputimes

  Ntests=10;

  match test_type:
    case 1:
      # TEST 1: concave arc
      cc=np.array([0, 0]); r=0.25   # circle defs; center "cc" and radius "r"
      a=np.array([0.25, 0])                        # first point of the circle
      b=np.array([0, 0.25])                        # last point of the circle
      v=np.array([[0.5, 0], [0.5, 0.5], [0, 0.5]]) # polygon vertices
      conv=0                                       # 0: concave arc

    case 2:
      # TEST 2: convex arc
      cc=np.array([.25, .25]); r=.25 # circle defs; center "cc" and radius "r"
      a=np.array([0.25, 0])                        # first point of the circle
      b=np.array([0, 0.25])                        # last point of the circle
      # polygon vertices
      v=np.array([[.4, .05], [.5, .25], [.45, .45],\
                  [.3, .5], [.1, .45]])
      conv=1                                       # 1: convex arc

    case 3:
      # TEST 3: concave arc
      cc=np.array([0, 0]); r=0.25   # circle defs; center "cc" and radius "r"
      a=np.array([0.25, 0])                        # first point of the circle
      b=np.array([0, 0.25])                        # last point of the circle
      v=np.array([[0.25, 0.2], [0.2, 0.25]])       # polygon vertices
      conv=0                                       # 0: concave arc

    case 4:
      # TEST 4 convex arc
      cc=np.array([0, 0]); r=0.25   # circle defs; center "cc" and radius "r"
      a=np.array([0.25, 0])                        # first point of the circle
      b=np.array([0, 0.25])                        # last point of the circle
      # polygon vertices
      v=np.array([[.4, .05], [.5, .25], [.45, .45], [.3, .5], [.1, .45]])
      conv=0                                       # 0: concave arc

# ...................... Compute cubature rules  ......................

  cpusC = []
  for k in range(Ntests):
    t = time.time()
    xyw = polygcirc(n,v,a,b,cc,r,conv,pos,False)
    cpusC += [time.time() - t]
  cpusM = [statistics.median(cpusC)]

  cpusC = []
  for k in range(Ntests):
    t = time.time()
    xyw,xywc,P,L,subs, momerr = polygcirc(n,v,a,b,cc,r,conv,pos)
    cpusC += [time.time() - t]
  cpusM += [statistics.median(cpusC)]

# .......................... Statistics  ...............................

  print('\n \t EXAMPLE         : %d' %test_type)
  print('\t ALG. DEG. PREC. : %d' %n)

  match pos:
    case 1:
      print('\t COMPRESSION     : LHDM')
    case _:
      print('\t COMPRESSION     : lsqnonneg')

  print('\t NODES FULL RULE : %d' %xyw[:,0].size)
  print('\t NODES COMP RULE : %d' %xywc[:,0].size)
  w=xyw[:,2]; iwneg=np.where(w < 0)[0]; L=len(iwneg)
  print('\t NEG WEIGHTS FULL: %d' %L)
  wc=xywc[:,2]; iwnegc=np.where(wc < 0)[0]; Lc=len(iwnegc);
  print('\t FULL RULE CPU   : %1.3e' %cpusM[0])
  print('\t FULL + COMP CPU : %1.3e' %cpusM[1])
  print('\t MOMENTS\' ERROR  : %1.3e' %momerr)

#..................... Nodes to copy in Matlab .......................

  # Cubature Full Nodes
#  print('[')
#  for i in xyw:
#    print('%.15f  %.15f %.15f;'  %(i[0], i[1], i[2]))
#  print(']')

  # Cubature Comprex Nodes
#  print('[')
#  for i in xywc:
#    print('%.15f  %.15f %.15f;'  %(i[0], i[1], i[2]))
#  print(']')

#  ......................... Plots  ..................................

  # ... plot parameters ...
  size_sides=2
  size_pt=4
  size_cmp=8
  size_in_sides=1

  # ... plot figure ...
  fig, axs = plt.subplots(figsize=(10, 10))
  vertices=np.vstack([a, v, b])
  tha = np.arctan2(a[1]-cc[1],a[0]-cc[0])
  thb = np.arctan2(b[1]-cc[1],b[0]-cc[0])
  LL=int(subs[:,0].size/4)

  # Plot circular quadrangles (divisive dotted lines)
  for kk in range(1,LL+1):
    CC = cc
    Li = 4*(kk-1)+1
    A=subs[Li-1]
    B=subs[Li]
    C=subs[Li+1]
    D=subs[Li+2]
    axs.plot([A[0], C[0]],[A[1], C[1]],'k--')
    axs.plot([B[0], D[0]],[B[1], D[1]],'k--')

  axs.plot(vertices[:,0], vertices[:,1], 'k-')  # edges
  Lin = 1
  L = [L]
  P = np.array(P)
  for ii in range(len(L)):
    Lfin=L[ii]+Lin-1
    xywr = xyw[Lin:Lfin]
    Lin = Lfin+1
    iid = (ii+1)%7

    match iid:
      case 1:
        axs.plot(xywr[:,1],xywr[:,2],'b+')
      case 2:
        axs.plot(xywr[:,1],xywr[:,2],'r+')
      case 3:
        axs.plot(xywr[:,1],xywr[:,2],'c+',markeredgecolor='c',\
                 markerfacecolor='c',markersize=size_pt)
      case 4:
        axs.plot(xywr[:,1],xywr[:,2],'m+')
      case 5:
        axs.plot(xywr[:,1],xywr[:,2],'g+')
      case 6:
        axs.plot(xywr[:,1],xywr[:,2],'y+')
      case _:
        axs.plot(xywr[:,1],xywr[:,2],'k+')

  if np.size(P,0)>0:
    axs.plot(P[:,0],P[:,1],'k--')
  axs.plot(xyw[:,0],xyw[:,1],'go',markeredgecolor='k', markersize=10)  #dots

  if conv == 1:
    if tha<thb:
      tha = tha+2*np.pi
    th = np.linspace(thb,tha,100)
  else:
    th = np.linspace(tha,thb,100)

  xP=cc[0]+r*np.cos(th)
  yP=cc[1]+r*np.sin(th)
  axs.plot(xP,yP,'k-')

#.......................Plots with sorted weights .....................

#  fig, wpl = plt.subplots(figsize=(10, 10))
#  wpl.plot(list(range(1,len(xyw[:,2])+1)), np.sort(xyw[:,2]),'o',\
#           color='red', markersize=7, markeredgecolor='midnightblue')
#
#  fig, wplc = plt.subplots(figsize=(10, 10))
#  wplc.plot(list(range(1,len(xywc[:,2])+1)), np.sort(xywc[:,2]),'o',\
#            color='red', markersize=7, markeredgecolor='midnightblue')

  return
\end{minted}