function [pts_wam,quadrangles_storage]=wam_polygon_quad(deg,vertices)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% The following routine computes a WAM over a polygon with vertices
% stored in "vertices", of degree "deg".
%
% The routine consists in a quadrangulation of the polygon and then
% determine a WAM in each quadrangle.
%
% The union of these WAM is a WAM in the polygon.
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% deg: WAM degree.
% vertices : abscissae of the vertices of the polygonal region;
%            it consists of a L x 2 matrix of points describing the
%           vertices.
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% wam_pts: wam points at degree "deg", it is a matrix "N x 2", where
%      "wam_pts(k,:)" is a vector containing the coordinates of the k-th
%      point of the wam.
%
% quadrangles_storage: cell whose k-th component are the vertices of the
%      k-th quadrangle (or triangle in case of degenerate quadrangles).
%--------------------------------------------------------------------------
%% Copyright (C) 2009- Mariano Gentile, Alvise Sommariva, Marco Vianello.
%%
%% 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
%%
%% Author:
%%          Mariano Gentile
%%          Alvise Sommariva <alvise@math.unipd.it>
%%          Marco Vianello   <marcov@math.unipd.it>
%%
%% Date    : -2009
%% Modified: January 3, 2026
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
% Subroutines required
%--------------------------------------------------------------------------
% 1. quadrangulation_algorithm
% 2. wam_quadrangles
%--------------------------------------------------------------------------

[vertices_sequence,b]=quadrangulation_algorithm(vertices);

pts_wam=[];
v=[];

L1=1; L2=0;

quadrangles_storage={};

for index=1:length(b)

    L2=b(index)+L2;
    quadrangleL=[vertices_sequence(L1:L2,:); vertices_sequence(L1,:)];
    L1=L2+1;
    pts_loc=wam_quadrangles(quadrangleL,deg);
    pts_wam=[pts_wam; pts_loc];
    
    quadrangles_storage{end+1}=quadrangleL;

end











%--------------------------------------------------------------------------
% COMPUTING WAM.
%--------------------------------------------------------------------------

function pts=wam_quadrangles(polygon_vertices,deg)

%--------------------------------------------------------------------------
% OBJECT.
%--------------------------------------------------------------------------
% COMPUTE MINIMAL TRIANGULATION/QUADRANGULATION ON CONVEX DOMAINS AND A WAM
% OF DEGREE deg ON SUCH POLYGON.
%--------------------------------------------------------------------------
% INPUT.
%--------------------------------------------------------------------------
% polygon_vertices: N x 2 matrix of cartesian coordinates of vertices of
%       the polygon without the assumption that first and last row must be
%        equal;
%
% deg:  wam degree.
%--------------------------------------------------------------------------
% OUTPUT.
%--------------------------------------------------------------------------
% pts: WAM points stored as M x 2 matrix.
%--------------------------------------------------------------------------
% SUBROUTINES.
%--------------------------------------------------------------------------
% 1. triangular_wam_2009 (external to this file).
% 2. triangle_map (external to this file).
% 3. pdpts (external to this file).
% 4. gausslegendre (external to this file).
% 5. tens_gausslobatto (external to this file).
%
% PS: Notice the only one routine between 3., 4. and 5. will be used,
% depending on user's preference.
%
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%% Copyright (C) 2007-2010 Len Bos, Alvise Sommariva, Marco Vianello.
%%
%% 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:
%%
%%          Len Bos          <leonardpeter.bos@univr.it>
%%          Alvise Sommariva <alvise@euler.math.unipd.it>
%%          Marco Vianello   <marcov@euler.math.unipd.it>
%%
%% Date: June 07, 2010.
%--------------------------------------------------------------------------

% CHECKING IF THE FIRST AND LAST VERTEX COINCIDE, OTHERWISE ADD A LAST
% VERTEX EQUAL TO THE FIRST ONE.
v1=polygon_vertices(1,:);
vend=polygon_vertices(end,:);

if (norm(v1-vend) > 0)
    polygon_vertices=[polygon_vertices; polygon_vertices(1,:)];
end

lx=size(polygon_vertices,1)-1; % number of total sides.

% CHECKING IF ONE POLYGON OF THE DECOMPOSITION WILL BE A TRIANGLE AND
% COMPUTE REFERENCE WAM.
if rem(lx,2) == 0
    pts_tr=[];
else
    pts_tr=triangular_wam_2009(deg);
end

% CHECKING IF ONE POLYGON OF THE DECOMPOSITION WILL BE A QUADRANGLE AND
% COMPUTE REFERENCE WAM.
if (lx >= 4)

    % ** USING PADUA POINTS. **
    % THE FACTOR 2 ON deg IS DUE TO HOW THE COONS MAPPINGS COMPUTE THE
    % WAM (BILINEAR MAP).
    % pts_sq = pdpts(2*deg,[-1,1,-1,1]);

    % ** USING TENSOR PRODUCT GAUSS-LOBATTO. **
    % THE FACTOR 1 ON deg IS DUE TO HOW THE COONS MAPPINGS COMPUTE THE
    % WAM (BILINEAR MAP).
    % FIXED deg, THE ROUTINE tens_gausslobatto PRODUCES A TENSORIAL RULE
    % BASED ON (deg+2) GAUSS-JACOBI-LOBATTO NODES ON THE INTERVAL [-1,1].
    % HENCE THE DIMENSION OF pts_sq IS (deg+2)^2.
    % THE POINTS ARE RELATIVE TO THE WEIGHT FUNCTION w(t)=(1-t)^a(1+t)^b.
    % BY SETTING tens_jaclob(deg,a,b).
    pts_sq=tens_gausslobatto(deg);
else
    pts_sq=[];
end

% COMPUTE THE NUMBER OF SUBPOLYGONS WHOSE UNION MAKES THE INITIAL POLYGON.
subdv=ceil(lx/2)-1;

% COMPUTING WAM BY UNION OF WAMS ON QUADRANGLES AND AT MOST A TRIANGLE.
pts=[];
iinit(1)=2; ifin(1)=min(4,lx);
for ii=1:subdv
    nsides=ifin(ii)-iinit(ii)+2;
    if (nsides == 4)
        ivv=[1; (iinit(ii):ifin(ii))'];
        vv=polygon_vertices(ivv,:);
        pts_vv=coons(vv,pts_sq);
    else
        % THE LAST QUADRANGLE OF THE DECOMPOSITION OF THE INITIAL
        % POLYGON THAT IS NOT A TRIANGLE IS ACTUALLY A TRIANGLE.
        ivv=[1; lx-1; lx];
        vv=polygon_vertices(ivv,:);
        pts_vv=triangle_map(vv,pts_tr);
    end
    pts=[pts; pts_vv];
    iinit(ii+1)=ifin(ii);
    ifin(ii+1)=min(ifin(ii)+2,lx);
end







function pts_out=triangle_map(triangle_vertices,pts_in)

%--------------------------------------------------------------------------
% OBJECT.
%--------------------------------------------------------------------------
% THIS FUNCTION MAPS POINTS pts_in OF THE UNIT TRIANGLE IN THE TRIANGLE
% WITH (ORDERED) VERTICES triangle_vertices, OBTAINING THE POINTS pts_out.
%--------------------------------------------------------------------------
% INPUTS.
%--------------------------------------------------------------------------
% triangle_vertices: vertices of a triangle.
%        It is not mandatory that first and last vertex coincide.
%
% pts_in: points on the reference triangle (0,0),(1,0),(0,1).
%--------------------------------------------------------------------------
% OUTPUT.
%--------------------------------------------------------------------------
% pts_out: the points "pts_in" are mapped on the triangle
%          "triangle_vertices" by a classical map.
%--------------------------------------------------------------------------

Q1=triangle_vertices(1,:);
Q2=triangle_vertices(2,:);
Q3=triangle_vertices(3,:);

b=Q1;
P=[(Q2-Q1)' (Q3-Q1)'];

pts_out_pre=(P*pts_in')';
pts_out(:,1)=b(1)+pts_out_pre(:,1);
pts_out(:,2)=b(2)+pts_out_pre(:,2);











function [vertices_sequence,b,A]=quadrangulation_algorithm(vertices)

%--------------------------------------------------------------------------
% Object.
%--------------------------------------------------------------------------
% This routine subdivides a normal polygon (described counterwise through
% its vertices vertices) into a sequence of quadrangles or triangles whose
% union is the initial polygon.
%
% Being the polygon normal, it must not have self-intersections.
%
% 1. the quadrangles/triangles are stored in the matrix "vertices_sequence"
% 2. the vector "b" describes if the subpolygon is a quadrangle or a
% triangle.
%
% The i-th subdomain has "b(i)" vertices that are the points of
% "vertices_sequence" with indices ranging from sum(b(1:I-1))+1 to
% sum(b(1:I-1))+b(I) with sum(b(1:0))=0.
% to clearify this representation, see the example below.
%
% Another way of storing the points is addressed in the matrix "A" (less
% handy, but proposed in gentile degree thesis).
%--------------------------------------------------------------------------
% Example.
%--------------------------------------------------------------------------
% Consider the domain
%
%     polygon_sides=[ 1 1; 6 2; 7 4; 10 3; 9 6; 6 7; 4 5; 1 6; 1 1];
%
% we get:
%
%     vertices_sequence = [
%         7.0000    4.0000
%        10.0000    3.0000
%         9.0000    6.0000
%         8.1429    6.2857
%         4.0000    5.0000
%         1.0000    6.0000
%         1.0000    2.0000
%         1.0000    2.0000
%         1.0000    1.0000
%         6.0000    2.0000
%         8.1429    6.2857
%         1.0000    2.0000
%         8.1429    6.2857
%         6.0000    7.0000
%         ];
%
%     b = [4     3     4     3];
%
% The first 4 rows of vertices_sequence represent the first quadrangle, the
% rows from 5 to 7 the first triangle, the rows from 8 to 11 the second
% quadrangle, the rows from 12 to 14 the second and last quadrangle.
%--------------------------------------------------------------------------
%% Copyright (C) 2010-2011
%% Mariano Gentile, Alvise Sommariva, Marco Vianello.
%%
%% 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
%%
%% Author:
%%          Mariano Gentile <gentimar@hotmail.it>
%%          Alvise Sommariva  <alvise@math.unipd.it>
%%          Marco Vianello    <marcov@math.unipd.it>
%%
%% Date       : December 26, 2010 - April 1, 2011.
%% Modified on: January 3, 2026.
%--------------------------------------------------------------------------

% LAST VERTEX MUST NOT BE REPEATED. ADJUSTING POSSIBLE PROBLEMS.
pts_init=vertices(1,:);
pts_end=vertices(end,:);

if norm(pts_init-pts_end) == 0
    vertices=vertices(1:end-1,:);
end

% QUADRANGULATION CODE STARTS HERE.

c=vertices;

A_nc=c;
b_nc(1)=size(c,1);
k_A_1=1; k_b_1=1;
t=1;
while t==1
    clear c_temp;
    for k=1:1:2
        for h=1:1:b_nc(1)
            c_temp(h,k)=A_nc(h,k);
        end
    end
    clear A_temp; clear b_temp;
    [A_temp,b_temp]=subdivision_concave_polygons(c_temp);
    % c_temp %%%
    %A_temp %%%
    %b_temp %%%
    if size(A_temp,2)==2
        for h=1:1:b_temp
            A_1(h,k_A_1)=A_temp(h,1);
            A_1(h,k_A_1+1)=A_temp(h,2);
        end
        %A_1 %%%
        k_A_1=k_A_1+2;
        b_1(k_b_1)=b_temp;
        k_b_1=k_b_1+1;
    end
    clear s_A_nc;
    clear s_A_temp;
    s_A_nc=size(A_nc,2);
    s_A_temp=size(A_temp,2);
    if (s_A_nc==2 && s_A_temp==2)
        t=0;
    end
    if (s_A_nc==2 && s_A_temp~=2)
        clear A_nc; clear b_nc;
        A_nc=A_temp;
        b_nc=b_temp;
    end
    if (s_A_nc~=2 && s_A_temp==2)
        clear A_nc_temp; clear b_nc_temp; clear k_b_nc_temp;
        A_nc_temp=A_nc;
        b_nc_temp=b_nc;
        clear A_nc; clear b_nc;
        k_b_nc_temp=1;
        for h=3:2:size(A_nc_temp,2)-1
            k_b_nc_temp=k_b_nc_temp+1;
            b_nc(k_b_nc_temp-1)=b_nc_temp(k_b_nc_temp);
            for k=1:1:b_nc_temp(k_b_nc_temp)
                A_nc(k,h-2)=A_nc_temp(k,h);
                A_nc(k,h-1)=A_nc_temp(k,h+1);
            end
        end
    end
    if (s_A_nc~=2 && s_A_temp~=2)
        clear A_nc_temp; clear b_nc_temp; clear k_b_nc_temp;
        A_nc_temp=A_nc;
        b_nc_temp=b_nc;
        clear A_nc; clear b_nc;
        k_b_nc_temp=1;
        for h=3:2:size(A_nc_temp,2)-1
            k_b_nc_temp=k_b_nc_temp+1;
            b_nc(k_b_nc_temp-1)=b_nc_temp(k_b_nc_temp);
            for k=1:1:b_nc_temp(k_b_nc_temp)
                A_nc(k,h-2)=A_nc_temp(k,h);
                A_nc(k,h-1)=A_nc_temp(k,h+1);
            end
        end
        f_A=size(A_nc_temp,2)-2;
        f_b=length(b_nc_temp)-1;
        for k=1:1:max(b_temp(1),b_temp(2))
            for h=1:1:4
                A_nc(k,h+f_A)=A_temp(k,h);
            end
        end
        b_nc(f_b+1)=b_temp(1);
        b_nc(f_b+2)=b_temp(2);
    end
end

k_A_1=1; k_A=0; k_b=0;
% A_1 %%%
% b_1 %%%
for k=1:1:length(b_1)
    clear c_convex;
    for l=1:1:b_1(k)
        c_convex(l,1)=A_1(l,k_A_1);
        c_convex(l,2)=A_1(l,k_A_1+1);
    end
    k_A_1=k_A_1+2;
    clear A_convex; clear b_convex;
    % c_convex %%%
    [A_convex,b_convex]=subdivision_convex_polygons(c_convex);
    for m=1:1:length(b_convex)
        b(k_b+m)=b_convex(m);
    end
    k_b=k_b+length(b_convex);

    for m=1:1:size(A_convex,2)
        for l=1:1:size(A_convex,1)
            A(l,k_A+m)=A_convex(l,m);
        end
    end
    k_A=k_A+size(A_convex,2);
end

s_A=size(A,2);   %%% eliminazione poligoni "senza area"
clear k_b;
k_b=0;
for k=1:2:s_A-1
    k_b=k_b+1;
    clear v_1; clear v_2; clear v_3; clear v_4;
    v_1(1)=A(1,k); v_1(2)=A(1,k+1); v_1(3)=0;
    v_2(1)=A(2,k); v_2(2)=A(2,k+1); v_2(3)=0;
    v_3(1)=A(3,k); v_3(2)=A(3,k+1); v_3(3)=0;
    v_4(1)=A(4,k); v_4(2)=A(4,k+1); v_4(3)=0;
    s_b=b(k_b);
    if s_b==3
        if norm(cross(v_3-v_2,v_2-v_1))==0
            b(k_b)=0;
        end
    end
    if s_b==4
        pv_1=norm(cross(v_1-v_4,v_2-v_1));
        pv_2=norm(cross(v_2-v_1,v_3-v_2));
        pv_3=norm(cross(v_4-v_3,v_3-v_2));
        pv_4=norm(cross(v_4-v_3,v_1-v_4));
        if pv_1==0
            b(k_b)=3;
            A(1,k)=A(2,k); A(1,k+1)=A(2,k+1);
            A(2,k)=A(3,k); A(2,k+1)=A(3,k+1);
            A(3,k)=A(4,k); A(3,k+1)=A(4,k+1);
            A(4,k)=0; A(4,k+1)=0;
        end
        if pv_2==0
            b(k_b)=3;
            A(2,k)=A(3,k); A(2,k+1)=A(3,k+1);
            A(3,k)=A(4,k); A(3,k+1)=A(4,k+1);
            A(4,k)=0; A(4,k+1)=0;
        end
        if pv_3==0
            b(k_b)=3;
            A(3,k)=A(4,k); A(3,k+1)=A(4,k+1);
            A(4,k)=0; A(4,k+1)=0;
        end
        if pv_4==0
            b(k_b)=3;
            A(4,k)=0; A(4,k+1)=0;
        end
        if pv_1+pv_2==0
            b(k_b)=0;
        end
    end
end
A_fin=A; s_A_fin=size(A_fin,2);
b_fin=b;
clear A; clear b;
k_b=1;
k_A=1;
k_b_fin=0;
for k=1:2:s_A_fin-1
    k_b_fin=k_b_fin+1;
    if b_fin(k_b_fin)~=0
        A(:,k_A)=A_fin(:,k);
        A(:,k_A+1)=A_fin(:,k+1);
        b(k_b)=b_fin(k_b_fin);
        k_A=k_A+2;
        k_b=k_b+1;
    end
end

vertices_sequence=polygon_decoder(A,b);












%--------------------------------------------------------------------------
% ATTACHED FUNCTIONS.
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
% 1. subdivision_convex_polygons
%--------------------------------------------------------------------------

function [A,b]=subdivision_convex_polygons(c)

% THIS ROUTINE SUBDIVIDES A CONVEX POLYGON (DESCRIBED COUNTERWISE THROUGH
% ITS VERTICES) INTO A SEQUENCE OF QUADRANGLES OR TRIANGLES WHOSE UNION IS
% THE INITIAL POLYGON.

% THE QUADRANGLES/TRIANGLES ARE STORED IN THE MATRIX "A" AND THE VECTOR "b"
% DESCRIBES IF THE SUBPOLYGON IS A QUADRANGLE OR A TRIANGLE.

%--------------------------------------------------------------------------
%% Copyright (C) 2010
%% Mariano Gentile, Alvise Sommariva, Marco Vianello.
%%
%% 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
%%
%% Author:
%%          Mariano Gentile
%%          Alvise Sommariva  <alvise@math.unipd.it>
%%          Marco Vianello    <marcov@math.unipd.it>
%%
%% Date: November 13, 2010.
%--------------------------------------------------------------------------


if size(c,2) ~= 2  % CHECKING POSSIBLE ERRORS IN THE POLYGON DESCRIPTION.
    print('error');
end

k=1; c_1=c; c_2=c_1; k_b=1;

while size(c_1,1) >= 4   % SUBDIVISION ALGORITHM.
    clear c_2;
    for l=1:1:4
        A(l,k)=c_1(l,1);
        A(l,k+1)=c_1(l,2);
    end
    b(k_b)=4;
    c_2(1,1)=c_1(1,1);
    c_2(1,2)=c_1(1,2);
    for m=4:1:size(c_1,1)
        c_2(m-2,1)=c_1(m,1);
        c_2(m-2,2)=c_1(m,2);
    end
    c_1=c_2;
    k=k+2;
    k_b=k_b+1;
end

if size(c_1,1) == 3   % POSSIBLE TRIANGLE.
    b(k_b)=3;
    for t=1:1:3
        A(t,k)=c_1(t,1);
        A(t,k+1)=c_1(t,2);
    end
end





%--------------------------------------------------------------------------
% 2. subdivision_concave_polygons.
%--------------------------------------------------------------------------

function [A,b]=subdivision_concave_polygons(c)

% divisione di poligoni concavi in quadrilateri e triangoli

n=size(c,1);
m=0;
k=1;
while (m==0 && k~=n+1)
    v_1(1)=c(mod_1(k-1,n),1);
    v_1(2)=c(mod_1(k-1,n),2);
    v_1(3)=0;
    v_2(1)=c(mod_1(k,n),1);
    v_2(2)=c(mod_1(k,n),2);
    v_2(3)=0;
    v_3(1)=c(mod_1(k+1,n),1);
    v_3(2)=c(mod_1(k+1,n),2);
    v_3(3)=0;
    v=cross(v_2-v_1,v_3-v_2);
    v_s=-sign(v(3));
    m=floor((v_s+1)/2);
    k=k+1;
end

if m==0
    A=c;
    b=size(c,1);
end
if m==1
    k=k-1;   % angolo concavo
    v_1(1)=c(mod_1(k-1,n),1);
    v_1(2)=c(mod_1(k-1,n),2);
    v_2(1)=c(mod_1(k,n),1);
    v_2(2)=c(mod_1(k,n),2);
    z=(v_2(2)-v_1(2))/(v_2(1)-v_1(1));
    if z==-Inf
        z=Inf;
    end
    s=k+1;
    c_2=[0,0,0,realmax,0,0];
    while mod_1(s,n)~=mod_1(k,n)
        v_1_bis(1)=c(mod_1(s,n),1);
        v_1_bis(2)=c(mod_1(s,n),2);
        v_2_bis(1)=c(mod_1(s+1,n),1);
        v_2_bis(2)=c(mod_1(s+1,n),2);
        controllo=0;
        z_1=(v_2_bis(2)-v_1_bis(2))/(v_2_bis(1)-v_1_bis(1));
        if z_1==-Inf
            z_1=Inf;
        end
        if z_1==z
            if v_1_bis(2)-v_2(2)==z*(v_1_bis(1)-v_2(1))
                t=(v_1_bis(1)-v_1(1))/(v_2(1)-v_1(1));
                if t>1
                    if t<c_2(4)
                        c_2(3)=s;
                        c_2(4)=t;
                        c_2(5)=1;
                        c_2(6)=1;
                    end
                end
            end
        end
        if z_1~=z
            if z==Inf
                if (v_1_bis(1)==v_2(1) || (v_1_bis(1)-v_2(1))*...
                        (v_2_bis(1)-v_2(1))<=0)
                    controllo=1;
                end
            end
            if z~=Inf
                if ((v_1_bis(2)-v_2(2)-(z*(v_1_bis(1)-v_2(1))))*...
                        (v_2_bis(2)-v_2(2)-(z*(v_2_bis(1)-v_2(1)))))<=0
                    controllo=1;
                end
            end
            if controllo==1
                W=[v_2_bis(1)-v_1_bis(1), v_1(1)-v_2(1); v_2_bis(2)-...
                    v_1_bis(2), v_1(2)-v_2(2)];
                V=[v_1(1)-v_1_bis(1); v_1(2)-v_1_bis(2)];
                X=W\V;
                tol=1e-10;
                if X(1)>tol && X(1)<1
                    if X(2)>1
                        if X(2)<c_2(4)
                            c_2(1)=v_1_bis(1)+X(1)*(v_2_bis(1)-v_1_bis(1));
                            c_2(2)=v_1_bis(2)+X(1)*(v_2_bis(2)-v_1_bis(2));
                            c_2(3)=s;
                            c_2(4)=X(2);
                            c_2(5)=0;
                            c_2(6)=0;
                        end
                    end
                end
                if X(1)<=tol && X(1)>=-tol
                    if X(2)>1
                        if X(2)<c_2(4)
                            c_2(3)=s;
                            c_2(4)=X(2);
                            c_2(5)=1;
                            c_2(6)=0;
                        end
                    end
                end
            end
        end
        s=s+1;
    end

    k_1=k;
    t=1;
    b(1)=0;
    while mod_1(k_1,n)~=mod_1(c_2(3),n)
        A(t,1)=c(mod_1(k_1,n),1);
        A(t,2)=c(mod_1(k_1,n),2);
        b(1)=b(1)+1;
        t=t+1;
        k_1=k_1+1;
    end
    if c_2(5)==0
        A(t,1)=c(mod_1(c_2(3),n),1);
        A(t,2)=c(mod_1(c_2(3),n),2);
        b(1)=b(1)+1;
        A(t+1,1)=c_2(1);
        A(t+1,2)=c_2(2);
        b(1)=b(1)+1;
    end
    if c_2(5)==1
        if z==Inf
            if c(mod_1(k,n),1)~=c(mod_1(c_2(3)-1,n),1)
                A(t,1)=c(mod_1(c_2(3),n),1);
                A(t,2)=c(mod_1(c_2(3),n),2);
                b(1)=b(1)+1;
            end
        end
        if z~=Inf
            if c(mod_1(c_2(3),n),2)~=((c(mod_1(c_2(3)-1,n),2)-...
                    c(mod_1(k,n),2))/(c(mod_1(c_2(3)-1,n),1)-...
                    c(mod_1(k,n),1)))*(c(mod_1(c_2(3),n),1)-...
                    c(mod_1(k,n),1))+c(mod_1(k,n),2)
                A(t,1)=c(mod_1(c_2(3),n),1);
                A(t,2)=c(mod_1(c_2(3),n),2);
                b(1)=b(1)+1;
            end
        end
    end

    s_1=c_2(3);
    t=1;
    b(2)=0;
    if c_2(5)==0
        A(t,3)=c_2(1);
        A(t,4)=c_2(2);
        b(2)=b(2)+1;
        t=t+1;
        s_1=s_1+1;
    end
    if c_2(5)==1
        if z==Inf
            if c(mod_1(k,n),1)==c(mod_1(s_1,n),1) && ...
                    c(mod_1(k,n),1)==c(mod_1(s_1+1,n),1)
                s_1=s_1+1;
            end
        end
        if z~=Inf
            if c(mod_1(s_1,n),2)==((c(mod_1(s_1+1,n),2)- ...
                    c(mod_1(k,n),2))/(c(mod_1(s_1+1,n),1)- ...
                    c(mod_1(k,n),1)))*(c(mod_1(s_1,n),1)- ...
                    c(mod_1(k,n),1))+c(mod_1(k,n),2)
                s_1=s_1+1;
            end
        end
    end
    while mod_1(s_1,n)~=mod_1(k,n)
        A(t,3)=c(mod_1(s_1,n),1);
        A(t,4)=c(mod_1(s_1,n),2);
        b(2)=b(2)+1;
        t=t+1;
        s_1=s_1+1;
    end
end





%--------------------------------------------------------------------------
% 3. mod_1.
%--------------------------------------------------------------------------

function [a]=mod_1(b,c)

a=mod(b,c);

if a==0
    a=c;
end




%--------------------------------------------------------------------------
% 4. polygon_decoder.
%--------------------------------------------------------------------------

function vertices=polygon_decoder(A,b)


k_b=1;
ntri=0;
nquad=0;
v=[];

for k=1:2:(size(A,2)-1)
    clear xloc; clear yloc;
    for h=1:1:3
        xloc(h,1)=A(h,k);
        yloc(h,1)=A(h,k+1);
    end
    if b(k_b)==4
        xloc(4,1)=A(4,k);
        yloc(4,1)=A(4,k+1);
        nquad=nquad+1;
    else
        ntri=ntri+1;
    end

    vloc=[xloc yloc];
    v=[v; vloc];

    k_b=k_b+1;

end

vertices=v;




%
%
% function res=isleft(P1,P2,P3)
%
% %--------------------------------------------------------------------------
% % PURPOSE.
% %--------------------------------------------------------------------------
% % P1, P2, P3 ARE DIFFERENT POINTS OF THE PLANE. P3 CAN BE A SET OF "N"
% % VECTORS.
% % P1 AND P2 ARE STORED AS ROW VECTORS. P3 IS A "N x 2" MATRIX.
% % THIS ROUTINE SAYS IF P3(k,:) IS ON THE LEFT SIDE W.R.T. THE STRAIGHT
% % LINE "r" GOING FROM P1 TO P2 (WITH OBVIOUS ORIENTATION).
% % IF res(K)=1 THEN P3(k,:) IS ON THE LEFT OF THE STRAIGHTLINE "r".
% % IF res(K)=0 THEN P3(k,:) IS ON THE RIGHT  OF THE STRAIGHTLINE "r".
% % IF res(K)=0.5 IS OVER THE LINE "r".
% %--------------------------------------------------------------------------
%
% % STAY ON THE LEFT IS EQUIVALENT THAT THE CROSS PRODUCT VECTOR STAYS
% % "OVER" THE xy PLANE. WE MODIFY THE VECTORS INVOLVED SINCE THE BUILT-IN
% % MATLAB FUNCTION REQUIRES A VECTOR AT LEAST OF DIMENSION 3.
%
% v1=P2-P1;
% v2=P3-repmat(P2,size(P3,1),1);
%
% v1=[v1 0]; % v1, v2 ARE ROW VECTORS.
% v2=[v2 zeros(size(v2,1),1)];
%
% %--------------------------------------------------------------------------
% % CROSS PRODUCT CAN BE COMPUTED USING cross THAT ALLOWS ARGUMENTS "A" AND
% % "B" BE MATRICES OF THE SAME DIMENSION AND DOING THE CROSS PRODUCT BETWEEN
% % COMPONENTS.
% %
% % EXAMPLE OF USING cross:
% %
% % >>A=[1;0;0];
% % >>B=[0;1;0];
% % >>cross(A,B)
% % ans =
% %
% %     0
% %     0
% %     1
% %
% % >>C=[1;1;1];
% % >>D=[2;3;1];
% % >>cross(C,D)
% %
% % ans =
% %
% %    -2
% %     1
% %     1
% %
% % >>cross([A C],[B D])
% %
% % ans =
% %
% %     0    -2
% %     0     1
% %     1     1
% %--------------------------------------------------------------------------
%
% % MODIFY v1, v2 (COLUMN VECTORS INSTEAD OF ROW VECTORS).
% v1=v1';
% v2=v2';
%
% % ADJUSTING v1 SO THAT IT HAS THE SAME DIMENSION OF v2.
% v1=repmat(v1,1,size(v2,2));
%
% % COMPUTING "MULTIPLE" CROSS PRODUCTS "cross(v1(:,k),v2(:,k))"
% cp=cross(v1,v2); % EACH COLUMN IS OF THE FORM [0;0;c].
% res=(cp(3,:))';
% res=(sign(res)+1)/2;




function pts=tens_gausslobatto(deg)

xw=lobatto_jacobi(deg,0,0);
x=xw(:,1); % deg+2 Gauss-Lobatto points

[X,Y]=meshgrid(x);

pts=[X(:) Y(:)];






function xw=lobatto_jacobi(N,a,b)

% LOBATTO_JACOBI  Gauss-Lobatto quadrature rule for Jacobi weight
% function.
%
%    xw=LOBATTO_JACOBI(n,a,b) generates the (n+2)-point Gauss-
%    Lobatto rule for the Jacobi weight function on [-1,1] with
%    parameters a and b. The n+2 nodes are stored in the first
%    column, the n+2 weights in the second column, of the (n+2)x2
%    array xw. The call xw=LOBATTO_JACOBI(n,a) is the same as xw=
%    LOBATTO_JACOBI(n,a,a) and xw=LOBATTO_JACOBI(n) the same as
%    xw=LOBATTO_JACOBI(n,0,0).
%
%    REFERENCE: W. Gautschi,``High-order Gauss-Lobatto formulae'',
%    Numer. Algorithms 25 (2000), 213-222.
%

if nargin<2, a=0; end; if nargin<3, b=a; end
ab=r_jacobi(N+2,a,b);
ab(N+2,1)=(a-b)/(2*N+a+b+2);
ab(N+2,2)=4*(N+a+1)*(N+b+1)*(N+a+b+1)/((2*N+a+b+1)* ...
    (2*N+a+b+2)^2);
xw=gauss(N+2,ab);




function xw=gauss(N,ab)

% GAUSS Gauss quadrature rule.
%
%    Given a weight function w encoded by the nx2 array ab 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.
%

N0=size(ab,1);
if N0<N
    error('input array ab too short'),
end
J=zeros(N);
for n=1:N, J(n,n)=ab(n,1); end
for n=2:N
    J(n,n-1)=sqrt(ab(n,2));
    J(n-1,n)=J(n,n-1);
end
[V,D]=eig(J);
[D,I]=sort(diag(D));
V=V(:,I);
xw=[D ab(1,2)*V(1,:)'.^2];









function ab=r_jacobi(N,a,b)

% R_JACOBI 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 nargin<2, a=0; end;  if nargin<3, b=a; end
if((N<=0)||(a<=-1)||(b<=-1)) 
    error('parameter(s) out of range')
end
nu=(b-a)/(a+b+2);
mu=2^(a+b+1)*gamma(a+1)*gamma(b+1)/gamma(a+b+2);
if N==1, ab=[nu mu]; return, end 
N=N-1; n=1:N; nab=2*n+a+b;
A=[nu (b^2-a^2)*ones(1,N)./(nab.*(nab+2))];
n=2:N; nab=nab(n);
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=[A' [mu; B1; B']];



function pts_out=coons(vertices,pts_in)

%--------------------------------------------------------------------------
% INPUTS.
%---------
% vertices : VERTICES OF A QUADRANGLE. IT NEEDS THAT THE FIRST AND LAST
%            VERTEX ARE THE SAME.
% pts_in: POINTS IN THE UNIT SQUARE [-1,1] x [-1,1]. THEY ARE STORED AS AN
%         N x 2 MATRIX.
%--------------------------------------------------------------------------
% OUTPUTS.
%---------
% pts_out: POINTS IN THE QUADRANGLE HAVING vertices AS VERTICES. THEY ARE
%          STORED AS AN N x 2 MATRIX AND OBTAINED FROM pts_in BY COONS'
%          LINEAR BLENDING.
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%% Copyright (C) 2010 Alvise Sommariva, Marco Vianello.
%%
%% 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:
%%
%%          Alvise Sommariva <alvise@euler.math.unipd.it>
%%          Marco Vianello   <marcov@euler.math.unipd.it>
%%
%% Date: June 08, 2010.
%--------------------------------------------------------------------------

% CHECKING IF THE FIRST AND LAST VERTEX COHINCIDE, OTHERWISE ADD A LAST
% VERTEX EQUAL TO THE FIRST ONE.
v1=vertices(1,:);
vend=vertices(end,:);

if (norm(v1-vend) > 0)
    vertices=[vertices; vertices(1,:)];
end

% APPLYING COONS LINEAR BLENDING THAT MAPS THE UNIT SQUARE INTO THE
% QUADRANGLE "Q" WITH VERTICES vertices SO TO MAP POINTS pts_in OF [0,1]^2 TO
% "Q".
u=pts_in(:,1);
v=pts_in(:,2);

A=repmat(vertices(1,:),length(u),1);
B=repmat(vertices(2,:),length(u),1);
C=repmat(vertices(3,:),length(u),1);
D=repmat(vertices(4,:),length(u),1);

term1=A.*repmat( (1-u).*(1-v),1,2 );
term2=B.*repmat( (1+u).*(1-v),1,2 );
term3=C.*repmat( (1+u).*(1+v),1,2 );
term4=D.*repmat( (1-u).*(1+v),1,2 );

pts_out=(1/4)*(term1+term2+term3+term4);









%--------------------------------------------------------------------------
% triangular_wam_2009.
%--------------------------------------------------------------------------

function [pts,boundary_points]=triangular_wam_2009(deg)

%--------------------------------------------------------------------------
% PURPOSE.
%-----------
%
% FOR A GIVEN DEGREE deg IT COMPUTES THE POINTS pts OF A TRIANGULAR WAM.
% BOUNDARY POINTS boundary_points ARE ALSO SPECIFIED.
%
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% INPUT.
%-----------
%
% deg: WAM DEGREE ON THE UNIT TRINAGLE (SIMPLEX).
%
%--------------------------------------------------------------------------
% OUTPUT.
%-----------
%
% pts: WAM pts.
% boundary_points: POINTS ON THE BOUNDARY OF THE TRIANGLE.
%
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%% Copyright (C) 2007-2010 Len Bos, Alvise Sommariva, Marco Vianello.
%%
%% 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:
%%
%%          Len Bos          <leonardpeter.bos@univr.it>
%%          Alvise Sommariva <alvise@euler.math.unipd.it>
%%          Marco Vianello   <marcov@euler.math.unipd.it>
%%
%% Date: May 17, 2010.
%--------------------------------------------------------------------------

% COMPUTING WAM ON THE UNIT DISK.
n1=2*deg;
n2=n1;
j1=(0:1:deg-1);
j2=(0:1:deg);
[rho,theta]=meshgrid(cos(j1*pi/n1),j2*pi/n2);

% MAPPING TO THE UNIT TRIANGLE.
B=[rho(:).*cos(theta(:)) rho(:).*sin(theta(:))];

% MESH POINTS ON THE TRIANGLE.
meshT=[B(:,1).^2 B(:,2).^2];
nodes_x=meshT(:,1);
nodes_y=meshT(:,2);
pts=[nodes_x nodes_y; 0 0];

% COMPUTING BOUNDARY POINTS.
i1=find(rho == 1);
i2=find(theta == pi/2);
i3=find(theta == 0);
boundary_points=union(i1,i2);
boundary_points=union(boundary_points,i3);
boundary_points=[boundary_points; size(pts,1)];
